diff --git a/ComputingTasks.cc b/ComputingTasks.cc
index c7a297880932eca46b24863186db1c38eed39a76..4791059204664dec6102860fb64f02a6f3f058d6 100644
--- a/ComputingTasks.cc
+++ b/ComputingTasks.cc
@@ -174,11 +174,13 @@ PriorPosteriorFunctionStatement::writeOutput(ostream &output, const string &base
 VarModelStatement::VarModelStatement(const SymbolList &symbol_list_arg,
                                      const OptionsList &options_list_arg,
                                      const string &name_arg,
-                                     const SymbolTable &symbol_table_arg) :
+                                     const SymbolTable &symbol_table_arg,
+                                     const DynamicModel &dynamic_model_arg) :
   symbol_list(symbol_list_arg),
   options_list(options_list_arg),
   name(name_arg),
-  symbol_table(symbol_table_arg)
+  symbol_table(symbol_table_arg),
+  dynamic_model(dynamic_model_arg)
 {
 }
 
@@ -214,6 +216,7 @@ VarModelStatement::writeOutput(ostream &output, const string &basename, bool min
   symbol_list.writeOutput("options_.var.var_list_", output);
   output << "M_.var." << name << " = options_.var;" << endl
          << "clear options_.var;" << endl;
+  dynamic_model.writeVarExpectationFunctions(output, name);
 }
 
 StochSimulStatement::StochSimulStatement(const SymbolList &symbol_list_arg,
diff --git a/ComputingTasks.hh b/ComputingTasks.hh
index bba0efbfb404d1374352dbc31449387c9232cee4..e76ba3aeabcf017a65a485acd403d8ebb07cdcb3 100644
--- a/ComputingTasks.hh
+++ b/ComputingTasks.hh
@@ -117,11 +117,13 @@ private:
   const OptionsList options_list;
   const string &name;
   const SymbolTable &symbol_table;
+  const DynamicModel &dynamic_model;
 public:
   VarModelStatement(const SymbolList &symbol_list_arg,
                     const OptionsList &options_list_arg,
                     const string &name_arg,
-                    const SymbolTable &symbol_table_arg);
+                    const SymbolTable &symbol_table_arg,
+                    const DynamicModel &dynamic_model_arg);
   virtual void checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings);
   virtual void writeOutput(ostream &output, const string &basename, bool minimal_workspace) const;
 };
diff --git a/DynamicModel.cc b/DynamicModel.cc
index ce3ea1134b03a0e52b783dc25a496f64e9eb3ee2..bb9d0ba9ef493d79e5a039e1722df5a65ba37ca8 100644
--- a/DynamicModel.cc
+++ b/DynamicModel.cc
@@ -3141,6 +3141,14 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
   output << "];" << endl;
 }
 
+void
+DynamicModel::writeVarExpectationFunctions(ostream &output, const string &var_model_name) const
+{
+  // Write Var Forecast files
+  for (int eq = 0; eq < (int) equations.size(); eq++)
+    equations[eq]->writeVarExpectationFunction(output, var_model_name);
+}
+
 map<pair<int, pair<int, int > >, expr_t>
 DynamicModel::collect_first_order_derivatives_endogenous()
 {
diff --git a/DynamicModel.hh b/DynamicModel.hh
index 7f0177fbd0593e5bc8036338b145af75adcb9fb3..9f0f2994b0ccd0097b05d1f33b10e737b25f6efb 100644
--- a/DynamicModel.hh
+++ b/DynamicModel.hh
@@ -236,6 +236,9 @@ public:
   //! Writes model initialization and lead/lag incidence matrix to output
   void writeOutput(ostream &output, const string &basename, bool block, bool byte_code, bool use_dll, int order, bool estimation_present, bool compute_xrefs, bool julia) const;
 
+  //! Writes calls to functions that write functions that forecast VARs
+  void writeVarExpectationFunctions(ostream &output, const string &var_model_name) const;
+
   //! Return true if the hessian is equal to zero
   inline bool checkHessianZero() const;
 
diff --git a/DynareBison.yy b/DynareBison.yy
index 0b56a4a30eb0e14798e23fdf37eac3fa28f85272..55c058554fabb1d444d734e970e25de575fc3758 100644
--- a/DynareBison.yy
+++ b/DynareBison.yy
@@ -771,7 +771,7 @@ hand_side : '(' hand_side ')'
             { $$ = driver.add_power($1, $3); }
           | EXPECTATION '(' signed_integer ')''(' hand_side ')'
 	    { $$ = driver.add_expectation($3, $6); }
-          | VAR_EXPECTATION '(' hand_side COMMA MODEL_NAME EQUAL NAME ')'
+          | VAR_EXPECTATION '(' symbol COMMA MODEL_NAME EQUAL NAME ')'
             { $$ = driver.add_var_expectation($3, $7); }
           | MINUS hand_side %prec UMINUS
             { $$ = driver.add_uminus($2); }
diff --git a/ExprNode.cc b/ExprNode.cc
index fc768c5e38122ef9bd5f235952f2d0a7eba55b1d..97be9eb183af142f76f5193d7ef8b8575d3a75c0 100644
--- a/ExprNode.cc
+++ b/ExprNode.cc
@@ -194,6 +194,12 @@ ExprNode::compileExternalFunctionOutput(ostream &CompileCode, unsigned int &inst
   // Nothing to do
 }
 
+void
+ExprNode::writeVarExpectationFunction(ostream &output, const string &var_model_name) const
+{
+  // Nothing to do
+}
+
 VariableNode *
 ExprNode::createEndoLeadAuxiliaryVarForMyself(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
@@ -349,6 +355,11 @@ NumConstNode::collectDynamicVariables(SymbolType type_arg, set<pair<int, int> >
 {
 }
 
+void
+NumConstNode::writeVarExpectationFunction(ostream &output, const string &var_model_name) const
+{
+}
+
 pair<int, expr_t >
 NumConstNode::normalizeEquation(int var_endo, vector<pair<int, pair<expr_t, expr_t> > > &List_of_Op_RHS) const
 {
@@ -1394,6 +1405,11 @@ VariableNode::isNumConstNodeEqualTo(double value) const
   return false;
 }
 
+void
+VariableNode::writeVarExpectationFunction(ostream &output, const string &var_model_name) const
+{
+}
+
 bool
 VariableNode::isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const
 {
@@ -1668,6 +1684,17 @@ UnaryOpNode::composeDerivatives(expr_t darg, int deriv_id)
   exit(EXIT_FAILURE);
 }
 
+void
+UnaryOpNode::writeVarExpectationFunction(ostream &output, const string &var_model_name) const
+{
+  if (op_code == oVarExpectation && var_model_name == var_expectation_model_name)
+    {
+      VariableNode *varg = dynamic_cast<VariableNode *>(arg);
+      output << "writeVarExpectationFunction('" << var_expectation_model_name << "', '"
+             << datatree.symbol_table.getName(varg->symb_id) << "');" << endl;
+    }
+}
+
 expr_t
 UnaryOpNode::computeDerivative(int deriv_id)
 {
@@ -2017,15 +2044,13 @@ UnaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
       break;
     case oVarExpectation:
       VariableNode *varg = dynamic_cast<VariableNode *>(arg);
-      assert(varg != NULL);
-      assert(datatree.symbol_table.getType(varg->symb_id) == eEndogenous);
       if (IS_LATEX(output_type))
         output << "VAR" << LEFT_PAR(output_type)
                << datatree.symbol_table.getTeXName(varg->symb_id)
                << "_{t+1}" << RIGHT_PAR(output_type);
       else
-        output << "var_forecast('" << var_expectation_model_name << "', 1, y, '"
-               << datatree.symbol_table.getName(varg->symb_id) << "')";
+        output << "var_forecast_" << var_expectation_model_name << "_"
+               << datatree.symbol_table.getName(varg->symb_id) << "(y)";
       return;
     }
 
@@ -2385,6 +2410,9 @@ UnaryOpNode::buildSimilarUnaryOpNode(expr_t alt_arg, DataTree &alt_datatree) con
 expr_t
 UnaryOpNode::toStatic(DataTree &static_datatree) const
 {
+  if (op_code == oVarExpectation)
+    return static_datatree.AddVariable((dynamic_cast<VariableNode *>(arg))->symb_id);
+
   expr_t sarg = arg->toStatic(static_datatree);
   return buildSimilarUnaryOpNode(sarg, static_datatree);
 }
@@ -3095,6 +3123,13 @@ BinaryOpNode::containsExternalFunction() const
     || arg2->containsExternalFunction();
 }
 
+void
+BinaryOpNode::writeVarExpectationFunction(ostream &output, const string &var_model_name) const
+{
+  arg1->writeVarExpectationFunction(output, var_model_name);
+  arg2->writeVarExpectationFunction(output, var_model_name);
+}
+
 void
 BinaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
                           const temporary_terms_t &temporary_terms,
@@ -4234,6 +4269,14 @@ TrinaryOpNode::containsExternalFunction() const
     || arg3->containsExternalFunction();
 }
 
+void
+TrinaryOpNode::writeVarExpectationFunction(ostream &output, const string &var_model_name) const
+{
+  arg1->writeVarExpectationFunction(output, var_model_name);
+  arg2->writeVarExpectationFunction(output, var_model_name);
+  arg3->writeVarExpectationFunction(output, var_model_name);
+}
+
 void
 TrinaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
                            const temporary_terms_t &temporary_terms,
@@ -4899,6 +4942,14 @@ AbstractExternalFunctionNode::normalizeEquation(int var_endo, vector<pair<int, p
     return (make_pair(1, (expr_t) NULL));
 }
 
+void
+AbstractExternalFunctionNode::writeVarExpectationFunction(ostream &output, const string &var_model_name) const
+{
+  for (vector<expr_t>::const_iterator it = arguments.begin();
+       it != arguments.end(); it++)
+    (*it)->writeVarExpectationFunction(output, var_model_name);
+}
+
 void
 AbstractExternalFunctionNode::writeExternalFunctionArguments(ostream &output, ExprNodeOutputType output_type,
                                                              const temporary_terms_t &temporary_terms,
diff --git a/ExprNode.hh b/ExprNode.hh
index 8f6d356092eeae52c56335da3777fd9d1b01c4c9..19a2494b4e361c76a1daadc128cb28e25847f76f 100644
--- a/ExprNode.hh
+++ b/ExprNode.hh
@@ -212,6 +212,9 @@ public:
   //! returns true if the expr node contains an external function
   virtual bool containsExternalFunction() const = 0;
 
+  //! Writes the matlab/C function for a Var Expectation Node
+  virtual void writeVarExpectationFunction(ostream &output, const string &var_model_name) const = 0;
+
   //! Writes output of node (with no temporary terms and with "outside model" output type)
   void writeOutput(ostream &output) const;
 
@@ -501,6 +504,7 @@ public:
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual bool isNumConstNodeEqualTo(double value) const;
+  virtual void writeVarExpectationFunction(ostream &output, const string &var_model_name) const;
   virtual bool containsEndogenous(void) const;
   virtual bool containsExogenous() const;
   virtual bool isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const;
@@ -567,6 +571,7 @@ public:
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual bool isNumConstNodeEqualTo(double value) const;
+  virtual void writeVarExpectationFunction(ostream &output, const string &var_model_name) const;
   virtual bool containsEndogenous(void) const;
   virtual bool containsExogenous() const;
   virtual bool isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const;
@@ -655,6 +660,7 @@ public:
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual bool isNumConstNodeEqualTo(double value) const;
+  virtual void writeVarExpectationFunction(ostream &output, const string &var_model_name) const;
   virtual bool containsEndogenous(void) const;
   virtual bool containsExogenous() const;
   virtual bool isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const;
@@ -754,6 +760,7 @@ public:
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual bool isNumConstNodeEqualTo(double value) const;
+  virtual void writeVarExpectationFunction(ostream &output, const string &var_model_name) const;
   virtual bool containsEndogenous(void) const;
   virtual bool containsExogenous() const;
   virtual bool isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const;
@@ -835,6 +842,7 @@ public:
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual bool isNumConstNodeEqualTo(double value) const;
+  virtual void writeVarExpectationFunction(ostream &output, const string &var_model_name) const;
   virtual bool containsEndogenous(void) const;
   virtual bool containsExogenous() const;
   virtual bool isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const;
@@ -916,6 +924,7 @@ public:
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual bool isNumConstNodeEqualTo(double value) const;
+  virtual void writeVarExpectationFunction(ostream &output, const string &var_model_name) const;
   virtual bool containsEndogenous(void) const;
   virtual bool containsExogenous() const;
   virtual bool isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const;
diff --git a/ModFile.cc b/ModFile.cc
index 5719d0197fe9d8d27a65846b3a285f354c2be378..975b0d8af06b0d43fc8d05b6319db691a3eea650 100644
--- a/ModFile.cc
+++ b/ModFile.cc
@@ -809,9 +809,20 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool clear_glo
     }
 
   // Print statements
+  // Print var_model statements first so that the forecast functions can be used by following statements
+  //  for (vector<Statement *>::const_iterator it = statements.begin();
+  //       it != statements.end(); it++)
+  //    if (dynamic_cast<VarModelStatement *>(*it) != NULL)
+  //      (*it)->writeOutput(mOutputFile, basename, minimal_workspace);
+  //  dynamic_model.writeVarExpectationFunctions(mOutputFile);
+
   for (vector<Statement *>::const_iterator it = statements.begin();
        it != statements.end(); it++)
     {
+      // Don't print var_model statements again as they were already printed
+      //      if (dynamic_cast<VarModelStatement *>(*it) != NULL)
+      //        continue;
+
       (*it)->writeOutput(mOutputFile, basename, minimal_workspace);
 
       /* Special treatment for initval block: insert initial values for the
diff --git a/ParsingDriver.cc b/ParsingDriver.cc
index 6e03325c9666153eab6d515bf8cebef30a3bdeb1..250c328bb7dc52ea7802f86d2e5e7729cbaf4f85 100644
--- a/ParsingDriver.cc
+++ b/ParsingDriver.cc
@@ -1286,7 +1286,7 @@ ParsingDriver::var_model()
   if (it == options_list.string_options.end())
     error("You must pass the model_name option to the var_model statement.");
   const string *name = new string(it->second);
-  mod_file->addStatement(new VarModelStatement(symbol_list, options_list, *name, mod_file->symbol_table));
+  mod_file->addStatement(new VarModelStatement(symbol_list, options_list, *name, mod_file->symbol_table, mod_file->dynamic_model));
   symbol_list.clear();
   options_list.clear();
 }
@@ -1628,6 +1628,15 @@ ParsingDriver::copy_options(string *to_declaration_type, string *to_name1, strin
   delete from_subsample_name;
 }
 
+void
+ParsingDriver::check_symbol_is_endogenous(string *name)
+{
+  check_symbol_existence(*name);
+  int symb_id = mod_file->symbol_table.getID(*name);
+  if (mod_file->symbol_table.getType(symb_id) != eEndogenous)
+    error(*name + " is not endogenous");
+}
+
 void
 ParsingDriver::check_symbol_is_endogenous_or_exogenous(string *name)
 {
@@ -2387,9 +2396,10 @@ ParsingDriver::add_expectation(string *arg1, expr_t arg2)
 }
 
 expr_t
-ParsingDriver::add_var_expectation(expr_t arg1, string *arg2)
+ParsingDriver::add_var_expectation(string *arg1, string *arg2)
 {
-  expr_t varExpectationNode = data_tree->AddVarExpectation(arg1, *arg2);
+  check_symbol_is_endogenous(arg1);
+  expr_t varExpectationNode = data_tree->AddVarExpectation(add_model_variable(arg1), *arg2);
   return varExpectationNode;
 }
 
diff --git a/ParsingDriver.hh b/ParsingDriver.hh
index 9d5736db2bfb4257675e28d674edb1aeec8bd5bc..a66d2d3332d16890a9fcd69da3b8ebaa894a4fe9 100644
--- a/ParsingDriver.hh
+++ b/ParsingDriver.hh
@@ -87,6 +87,9 @@ private:
   //! Checks that a given symbol exists and is a parameter, and stops with an error message if it isn't
   void check_symbol_is_parameter(string *name);
 
+  //! Checks that a given symbol exists and is endogenous, and stops with an error message if it isn't
+  void check_symbol_is_endogenous(string *name);
+
   //! Checks that a given symbol was assigned within a Statement
   void check_symbol_is_statement_variable(string *name);
 
@@ -639,7 +642,7 @@ public:
   //! Writes token "E(arg1)(arg2)" to model tree
   expr_t add_expectation(string *arg1,  expr_t arg2);
   //! Writes token "VAR_EXPECTATION(arg1,arg2)" to model tree
-  expr_t add_var_expectation(expr_t arg1,  string *arg2);
+  expr_t add_var_expectation(string *arg1,  string *arg2);
   //! Writes token "exp(arg1)" to model tree
   expr_t add_exp(expr_t arg1);
   //! Writes token "log(arg1)" to model tree