diff --git a/src/ComputingTasks.cc b/src/ComputingTasks.cc
index c54f5cedd05cd909911fe691cdc330a359406eb0..7277a0c442b3f376988c308665e6e80e94c1aec6 100644
--- a/src/ComputingTasks.cc
+++ b/src/ComputingTasks.cc
@@ -528,30 +528,6 @@ VarModelStatement::writeJsonOutput(ostream &output) const
          << "\"model_name\": \"" << name << "\"}";
 }
 
-void
-VarModelStatement::createVarModelMFunction(ostream &output, const map<string, set<int>> &var_expectation_functions_to_write) const
-{
-  if (var_expectation_functions_to_write.find(name) == var_expectation_functions_to_write.end())
-    return;
-
-  stringstream ss;
-  set<int> horizons = var_expectation_functions_to_write.find(name)->second;
-  for (auto it = horizons.begin(); it != horizons.end(); it++)
-    {
-      if (it != horizons.begin())
-        ss << " ";
-      ss << *it;
-    }
-
-  output << "writeVarExpectationFunction('" << name << "', ";
-  if (horizons.size() > 1)
-    output << "[";
-  output << ss.rdbuf();
-  if (horizons.size() > 1)
-    output << "]";
-  output << ");" << endl;
-}
-
 VarEstimationStatement::VarEstimationStatement(OptionsList options_list_arg) :
   options_list(move(options_list_arg))
 {
@@ -5017,3 +4993,45 @@ GenerateIRFsStatement::writeJsonOutput(ostream &output) const
     }
   output << "}";
 }
+
+VarExpectationModelStatement::VarExpectationModelStatement(string model_name_arg, string variable_arg, string var_model_name_arg,
+                                                           string horizon_arg, expr_t discount_arg, const SymbolTable &symbol_table_arg) :
+  model_name{move(model_name_arg)}, variable{move(variable_arg)},
+  var_model_name{move(var_model_name_arg)}, horizon{move(horizon_arg)},
+  discount{discount_arg}, symbol_table{symbol_table_arg}
+{
+}
+
+void
+VarExpectationModelStatement::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const
+{
+  string mstruct = "M_.var_expectation." + model_name;
+  output << mstruct << ".var_model_name = '" << var_model_name << "';" << endl
+         << mstruct << ".horizon = " << horizon << ';' << endl;
+  auto disc_var = dynamic_cast<const VariableNode *>(discount);
+  if (disc_var)
+    output << mstruct << ".discount_index = " << disc_var->get_symb_id() << ';' << endl;
+  else
+    {
+      output << mstruct << ".discount_value = ";
+      discount->writeOutput(output);
+      output << ';' << endl;
+    }
+  output << mstruct << ".param_indices = [ ";
+  for (int param_id : aux_params_ids)
+    output << symbol_table.getTypeSpecificID(param_id) << ' ';
+  output << "];" << endl;
+}
+
+void
+VarExpectationModelStatement::writeJsonOutput(ostream &output) const
+{
+  output << "{\"statementName\": \"var_expectation_model\","
+         << "\"model_name\": \"" << model_name << "\", "
+         << "\"variable\": \"" << variable << "\", "
+         << "\"var_model_name\": \"" << var_model_name << "\", "
+         << "\"horizon\": \"" << horizon << "\", "
+         << "\"discount\": \"";
+  discount->writeOutput(output);
+  output << "\"}";
+}
diff --git a/src/ComputingTasks.hh b/src/ComputingTasks.hh
index dda5b1981f9ad711b2becad304afaed027fff6b7..f6e7a0948ec369f10dea1114c103f9ae252f312d 100644
--- a/src/ComputingTasks.hh
+++ b/src/ComputingTasks.hh
@@ -145,10 +145,11 @@ public:
 
 class VarModelStatement : public Statement
 {
-private:
+public:
   const SymbolList symbol_list;
   const OptionsList options_list;
   const string name;
+private:
   const SymbolTable &symbol_table;
   vector<int> eqnumber, lhs, orig_diff_var;
   vector<set<pair<int, int>>> rhs_by_eq; // rhs by equation
@@ -170,7 +171,6 @@ public:
   void checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings) override;
   void writeOutput(ostream &output, const string &basename, bool minimal_workspace) const override;
   void writeJsonOutput(ostream &output) const override;
-  void createVarModelMFunction(ostream &output, const map<string, set<int>> &var_expectation_functions_to_write) const;
 };
 
 class VarRestrictionsStatement : public Statement
@@ -1209,4 +1209,20 @@ public:
   void writeJsonOutput(ostream &output) const override;
 };
 
+class VarExpectationModelStatement : public Statement
+{
+public:
+  const string model_name, variable, var_model_name, horizon;
+  const expr_t discount;
+  const SymbolTable &symbol_table;
+  // List of generated auxiliary param ids, in variable-major order
+  vector<int> aux_params_ids; // TODO: move this to some new VarModelTable object
+
+public:
+  VarExpectationModelStatement(string model_name_arg, string variable_arg, string var_model_name_arg,
+                               string horizon_arg, expr_t discount_arg, const SymbolTable &symbol_table_arg);
+  void writeOutput(ostream &output, const string &basename, bool minimal_workspace) const override;
+  void writeJsonOutput(ostream &output) const override;
+};
+
 #endif
diff --git a/src/DataTree.cc b/src/DataTree.cc
index da7f1d18f5d3e5aeac9fbd586c3d73f057bdf4c2..36f7c08be25afa28e1be544de74291e87875a189 100644
--- a/src/DataTree.cc
+++ b/src/DataTree.cc
@@ -505,15 +505,13 @@ DataTree::AddExpectation(int iArg1, expr_t iArg2)
 }
 
 expr_t
-DataTree::AddVarExpectation(const int symb_id, const int forecast_horizon, const string &model_name)
+DataTree::AddVarExpectation(const string &model_name)
 {
-  assert(symbol_table.getType(symb_id) == SymbolType::endogenous);
-
-  auto it = var_expectation_node_map.find({ model_name, symb_id, forecast_horizon });
+  auto it = var_expectation_node_map.find(model_name);
   if (it != var_expectation_node_map.end())
     return it->second;
 
-  return new VarExpectationNode(*this, symb_id, forecast_horizon, model_name);
+  return new VarExpectationNode(*this, model_name);
 }
 
 expr_t
diff --git a/src/DataTree.hh b/src/DataTree.hh
index 3d51ae8f40ef47595cb96a390009055c3a3f9335..45c90723599d596e526c575b60651487538f1759 100644
--- a/src/DataTree.hh
+++ b/src/DataTree.hh
@@ -82,7 +82,7 @@ protected:
   external_function_node_map_t external_function_node_map;
 
   // (model_name, symb_id, forecast_horizon) -> VarExpectationNode
-  using var_expectation_node_map_t = map<tuple<string, int, int>, VarExpectationNode *>;
+  using var_expectation_node_map_t = map<string, VarExpectationNode *>;
   var_expectation_node_map_t var_expectation_node_map;
 
   // model_name -> PacExpectationNode
@@ -242,8 +242,8 @@ public:
   expr_t AddSteadyStateParam2ndDeriv(expr_t iArg1, int param1_symb_id, int param2_symb_id);
   //! Adds "arg1=arg2" to model tree
   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 "var_expectation(model_name)" to model tree
+  expr_t AddVarExpectation(const string &model_name);
   //! Adds pac_expectation command to model tree
   expr_t AddPacExpectation(const string &model_name);
   //! Adds a model local variable with its value
diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index b44bc6ba34e108c5f1346ecd369b77ae241ff280..0c4fccbe11834b12ec69972cdcc198df8b14774c 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -1515,34 +1515,6 @@ DynamicModel::writeDynamicMFile(const string &basename) const
   writeDynamicModel(basename, false, false);
 }
 
-void
-DynamicModel::fillVarExpectationFunctionsToWrite()
-{
-  for (auto &it : var_expectation_node_map)
-    var_expectation_functions_to_write[get<0>(it.first)].insert(get<2>(it.first));
-}
-
-map<string, set<int>>
-DynamicModel::getVarExpectationFunctionsToWrite() const
-{
-  return var_expectation_functions_to_write;
-}
-
-void
-DynamicModel::writeVarExpectationCalls(ostream &output) const
-{
-  for (const auto & it : var_expectation_functions_to_write)
-    {
-      int i = 0;
-      output << "dynamic_var_forecast_" << it.first << " = "
-             << "var_forecast_" << it.first << "(y);" << endl;
-
-      for (auto it1 = it.second.begin(); it1 != it.second.end(); it1++)
-        output << "dynamic_var_forecast_" << it.first << "_" << *it1 << " = "
-               << "dynamic_var_forecast_" << it.first << "(" << ++i << ", :);" << endl;
-    }
-}
-
 void
 DynamicModel::writeDynamicJuliaFile(const string &basename) const
 {
@@ -3628,13 +3600,6 @@ DynamicModel::getVarLhsDiffAndInfo(vector<int> &eqnumber, vector<bool> &diff,
     }
 }
 
-void
-DynamicModel::setVarExpectationIndices(map<string, pair<SymbolList, int>> &var_model_info)
-{
-  for (auto & equation : equations)
-    equation->setVarExpectationIndex(var_model_info);
-}
-
 void
 DynamicModel::addEquationsForVar(map<string, pair<SymbolList, int>> &var_model_info)
 {
@@ -6220,3 +6185,10 @@ DynamicModel::writeJsonParamsDerivativesFile(ostream &output, bool writeDetails)
          << ", " << third_derivs1_output.str()
          << "}";
 }
+
+void
+DynamicModel::substituteVarExpectation(const map<string, expr_t> &subst_table)
+{
+  for (auto & equation : equations)
+    equation = dynamic_cast<BinaryOpNode *>(equation->substituteVarExpectation(subst_table));
+}
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index 1ca5ed4b91f57dfb1c5fe95d5a40327e7b80f0cb..ed50d97baa593d92e0b428b13248e6f22e8be32c 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -100,8 +100,6 @@ private:
   void writeDynamicMFile(const string &basename) const;
   //! Writes dynamic model file (Julia version)
   void writeDynamicJuliaFile(const string &dynamic_basename) const;
-  //! Write Var Expectation calls
-  void writeVarExpectationCalls(ostream &output) const;
   //! Writes dynamic model file (C version)
   /*! \todo add third derivatives handling */
   void writeDynamicCFile(const string &basename, const int order) const;
@@ -320,8 +318,6 @@ public:
   void getVarLhsDiffAndInfo(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);
   //! Get Pac equation parameter info
@@ -435,16 +431,13 @@ public:
   //! Substitutes diff operator
   void substituteDiff(StaticModel &static_model, ExprNode::subst_table_t &diff_subst_table);
 
+  //! Substitute VarExpectation operators
+  void substituteVarExpectation(const map<string, expr_t> &subst_table);
+
   //! Table to undiff LHS variables for pac vector z
   void getUndiffLHSForPac(vector<int> &lhs, vector<expr_t> &lhs_expr_t, vector<bool> &diff, vector<int> &orig_diff_var,
                           vector<int> &eqnumber, map<string, int> &undiff, ExprNode::subst_table_t &diff_subst_table);
 
-  //! Fill var_expectation_functions_to_write
-  void fillVarExpectationFunctionsToWrite();
-
-  //! Get var_expectation_functions_to_write
-  map<string, set<int>> getVarExpectationFunctionsToWrite() const;
-
   //! Transforms the model by replacing trend variables with a 1
   void removeTrendVariableFromEquations();
 
diff --git a/src/DynareBison.yy b/src/DynareBison.yy
index 380d0bd4573405420d85bcbd3ec94f3582083edc..2f5b6caab235ff49ee0faf4ed5ccebc815afd5ae 100644
--- a/src/DynareBison.yy
+++ b/src/DynareBison.yy
@@ -118,7 +118,7 @@ class ParsingDriver;
 %token TEX RAMSEY_MODEL RAMSEY_POLICY RAMSEY_CONSTRAINTS PLANNER_DISCOUNT DISCRETIONARY_POLICY DISCRETIONARY_TOL
 %token <string> TEX_NAME
 %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 VALUES VAR VAREXO VAREXO_DET VARIABLE VAROBS VAREXOBS PREDETERMINED_VARIABLES VAR_EXPECTATION VAR_EXPECTATION_MODEL 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 PAC_MODEL QOQ YOY AOA UNDIFF PAC_EXPECTATION
 %left EQUAL_EQUAL EXCLAMATION_EQUAL
@@ -175,8 +175,9 @@ class ParsingDriver;
 %type <string> filename symbol vec_of_vec_value vec_value_list date_expr number
 %type <string> vec_value_1 vec_value signed_inf signed_number_w_inf
 %type <string> range vec_value_w_inf vec_value_1_w_inf
-%type <string> integer_range signed_integer_range sub_sampling_options list_sub_sampling_option
-%type <pair<string,string>> named_var_elem subsamples_eq_opt calibration_range
+%type <string> integer_range signed_integer_range
+%type <string> sub_sampling_options list_sub_sampling_option
+%type <pair<string,string>> named_var_elem subsamples_eq_opt calibration_range integer_range_w_inf
 %type <vector<pair<string,string>>> named_var named_var_1
 %type <SymbolType> change_type_arg
 %type <vector<string>> vec_str vec_str_1
@@ -298,6 +299,7 @@ statement : parameters
           | gmm_estimation
           | smm_estimation
           | shock_groups
+          | var_expectation_model
           ;
 
 dsample : DSAMPLE INT_NUMBER ';'
@@ -386,6 +388,29 @@ pac_model_options : o_pac_name
                     { driver.pac_model_undiff($3, $5); }
                   ;
 
+var_expectation_model : VAR_EXPECTATION_MODEL '(' var_expectation_model_options_list ')' ';'
+                        { driver.var_expectation_model(); }
+                      ;
+
+var_expectation_model_options_list : var_expectation_model_option
+                                   | var_expectation_model_options_list COMMA var_expectation_model_option
+                                   ;
+
+
+var_expectation_model_option : VARIABLE EQUAL symbol
+                               { driver.option_str("variable", $3); }
+                             | VAR_MODEL_NAME EQUAL symbol
+                               { driver.option_str("var_model_name", $3); }
+                             | HORIZON EQUAL INT_NUMBER
+                               { driver.option_num("horizon", $3); }
+                             | HORIZON EQUAL integer_range_w_inf
+                               { driver.option_num("horizon", "[ " + $3.first + ' ' + $3.second + " ]"); }
+                             | MODEL_NAME EQUAL symbol
+                               { driver.option_str("model_name", $3); }
+                             | DISCOUNT EQUAL expression
+                               { driver.var_expectation_model_discount = $3; }
+                             ;
+
 restrictions : RESTRICTIONS '(' symbol ')' ';' { driver.begin_VAR_restrictions(); }
                restrictions_list END ';' { driver.end_VAR_restrictions($3); }
              ;
@@ -912,10 +937,8 @@ hand_side : '(' hand_side ')'
             { $$ = driver.add_power($1, $3); }
           | EXPECTATION '(' signed_integer ')''(' hand_side ')'
 	    { $$ = driver.add_expectation($3, $6); }
-          | VAR_EXPECTATION '(' symbol COMMA MODEL_NAME EQUAL NAME ')'
-            { $$ = driver.add_var_expectation($3, "1", $7); }
-          | VAR_EXPECTATION '(' symbol COMMA INT_NUMBER COMMA MODEL_NAME EQUAL NAME ')'
-            { $$ = driver.add_var_expectation($3, $5, $9); }
+          | VAR_EXPECTATION '(' symbol ')'
+            { $$ = driver.add_var_expectation($3); }
           | PAC_EXPECTATION '(' symbol ')'
             { $$ = driver.add_pac_expectation($3); }
           | MINUS hand_side %prec UMINUS
@@ -3650,6 +3673,12 @@ range : symbol ':' symbol
 integer_range : INT_NUMBER ':' INT_NUMBER
                 { $$ = $1 + ':' + $3; }
 
+integer_range_w_inf : INT_NUMBER ':' INT_NUMBER
+                      { $$ = make_pair($1, $3); }
+                    | INT_NUMBER ':' INF_CONSTANT
+                      { $$ = make_pair($1, "Inf"); }
+                    ;
+
 signed_integer_range : signed_integer ':' signed_integer
                        { $$ = $1 + ':' + $3; }
                      | MINUS '(' signed_integer ':' signed_integer ')'
diff --git a/src/DynareFlex.ll b/src/DynareFlex.ll
index 9d7248920cfebbc68be993cf7f92effb33467f8e..127b9090f8b0c68e7b3d1ecf8aa71f014720a820 100644
--- a/src/DynareFlex.ll
+++ b/src/DynareFlex.ll
@@ -141,6 +141,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
 <INITIAL>simul {BEGIN DYNARE_STATEMENT; return token::SIMUL;}
 <INITIAL>stoch_simul {BEGIN DYNARE_STATEMENT; return token::STOCH_SIMUL;}
 <INITIAL>var_model {BEGIN DYNARE_STATEMENT; return token::VAR_MODEL;}
+<INITIAL>var_expectation_model {BEGIN DYNARE_STATEMENT; return token::VAR_EXPECTATION_MODEL;}
 <INITIAL>pac_model {BEGIN DYNARE_STATEMENT; return token::PAC_MODEL;}
 <INITIAL>dsample {BEGIN DYNARE_STATEMENT; return token::DSAMPLE;}
 <INITIAL>Sigma_e {BEGIN DYNARE_STATEMENT; sigma_e = 1; return token::SIGMA_E;}
@@ -675,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
 <DYNARE_STATEMENT>emas_drop {return token::EMAS_DROP; }
 <DYNARE_STATEMENT>emas_tolf {return token::EMAS_TOLF; }
 <DYNARE_STATEMENT>emas_max_iter {return token::EMAS_MAX_ITER; }
+<DYNARE_STATEMENT>variable {return token::VARIABLE;}
 
 <DYNARE_STATEMENT>[\$][^$]*[\$] {
   strtok(yytext+1, "$");
diff --git a/src/ExprNode.cc b/src/ExprNode.cc
index f2a6f3e47b74d065c416a834b0ca4fc32d08eea8..8a09162b42152a7f5d14df10e41e047212c3368f 100644
--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
@@ -547,6 +547,12 @@ NumConstNode::substituteAdl() const
   return const_cast<NumConstNode *>(this);
 }
 
+expr_t
+NumConstNode::substituteVarExpectation(const map<string, expr_t> &subst_table) const
+{
+  return const_cast<NumConstNode *>(this);
+}
+
 void
 NumConstNode::findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const
 {
@@ -643,11 +649,6 @@ NumConstNode::isInStaticForm() const
   return true;
 }
 
-void
-NumConstNode::setVarExpectationIndex(map<string, pair<SymbolList, int>> &var_model_info)
-{
-}
-
 void
 NumConstNode::walkPacParameters(bool &pac_encountered, pair<int, int> &lhs, set<pair<int, pair<int, int>>> &ec_params_and_vars, set<pair<int, pair<int, int>>> &ar_params_and_vars) const
 {
@@ -1419,6 +1420,12 @@ VariableNode::substituteAdl() const
   return const_cast<VariableNode *>(this);
 }
 
+expr_t
+VariableNode::substituteVarExpectation(const map<string, expr_t> &subst_table) const
+{
+  return const_cast<VariableNode *>(this);
+}
+
 void
 VariableNode::findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const
 {
@@ -1793,11 +1800,6 @@ VariableNode::isInStaticForm() const
   return lag == 0;
 }
 
-void
-VariableNode::setVarExpectationIndex(map<string, pair<SymbolList, int>> &var_model_info)
-{
-}
-
 void
 VariableNode::walkPacParameters(bool &pac_encountered, pair<int, int> &lhs, set<pair<int, pair<int, int>>> &ec_params_and_vars, set<pair<int, pair<int, int>>> &ar_params_and_vars) const
 {
@@ -3009,6 +3011,13 @@ UnaryOpNode::substituteAdl() const
   return retval;
 }
 
+expr_t
+UnaryOpNode::substituteVarExpectation(const map<string, expr_t> &subst_table) const
+{
+  expr_t argsubst = arg->substituteVarExpectation(subst_table);
+  return buildSimilarUnaryOpNode(argsubst, datatree);
+}
+
 int
 UnaryOpNode::countDiffs() const
 {
@@ -3390,12 +3399,6 @@ UnaryOpNode::isInStaticForm() const
     return arg->isInStaticForm();
 }
 
-void
-UnaryOpNode::setVarExpectationIndex(map<string, pair<SymbolList, int>> &var_model_info)
-{
-  arg->setVarExpectationIndex(var_model_info);
-}
-
 void
 UnaryOpNode::walkPacParameters(bool &pac_encountered, pair<int, int> &lhs, set<pair<int, pair<int, int>>> &ec_params_and_vars, set<pair<int, pair<int, int>>> &ar_params_and_vars) const
 {
@@ -4870,6 +4873,15 @@ BinaryOpNode::substituteAdl() const
   return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
 }
 
+
+expr_t
+BinaryOpNode::substituteVarExpectation(const map<string, expr_t> &subst_table) const
+{
+  expr_t arg1subst = arg1->substituteVarExpectation(subst_table);
+  expr_t arg2subst = arg2->substituteVarExpectation(subst_table);
+  return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
+}
+
 void
 BinaryOpNode::findUnaryOpNodesForAuxVarCreation(DataTree &static_datatree, diff_table_t &nodes) const
 {
@@ -4991,13 +5003,6 @@ BinaryOpNode::isInStaticForm() const
   return arg1->isInStaticForm() && arg2->isInStaticForm();
 }
 
-void
-BinaryOpNode::setVarExpectationIndex(map<string, pair<SymbolList, int>> &var_model_info)
-{
-  arg1->setVarExpectationIndex(var_model_info);
-  arg2->setVarExpectationIndex(var_model_info);
-}
-
 void
 BinaryOpNode::walkPacParametersHelper(const expr_t arg1, const expr_t arg2,
                                       pair<int, int> &lhs,
@@ -5779,6 +5784,16 @@ TrinaryOpNode::substituteAdl() const
   return buildSimilarTrinaryOpNode(arg1subst, arg2subst, arg3subst, datatree);
 }
 
+
+expr_t
+TrinaryOpNode::substituteVarExpectation(const map<string, expr_t> &subst_table) const
+{
+  expr_t arg1subst = arg1->substituteVarExpectation(subst_table);
+  expr_t arg2subst = arg2->substituteVarExpectation(subst_table);
+  expr_t arg3subst = arg3->substituteVarExpectation(subst_table);
+  return buildSimilarTrinaryOpNode(arg1subst, arg2subst, arg3subst, datatree);
+}
+
 void
 TrinaryOpNode::findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const
 {
@@ -5901,14 +5916,6 @@ TrinaryOpNode::isInStaticForm() const
   return arg1->isInStaticForm() && arg2->isInStaticForm() && arg3->isInStaticForm();
 }
 
-void
-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);
-}
-
 void
 TrinaryOpNode::walkPacParameters(bool &pac_encountered, pair<int, int> &lhs, set<pair<int, pair<int, int>>> &ec_params_and_vars, set<pair<int, pair<int, int>>> &ar_params_and_vars) const
 {
@@ -6214,6 +6221,15 @@ AbstractExternalFunctionNode::substituteAdl() const
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
+expr_t
+AbstractExternalFunctionNode::substituteVarExpectation(const map<string, expr_t> &subst_table) const
+{
+  vector<expr_t> arguments_subst;
+  for (auto argument : arguments)
+    arguments_subst.push_back(argument->substituteVarExpectation(subst_table));
+  return buildSimilarExternalFunctionNode(arguments_subst, datatree);
+}
+
 void
 AbstractExternalFunctionNode::findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const
 {
@@ -6398,13 +6414,6 @@ AbstractExternalFunctionNode::isInStaticForm() const
   return true;
 }
 
-void
-AbstractExternalFunctionNode::setVarExpectationIndex(map<string, pair<SymbolList, int>> &var_model_info)
-{
-  for (auto argument : arguments)
-    argument->setVarExpectationIndex(var_model_info);
-}
-
 void
 AbstractExternalFunctionNode::walkPacParameters(bool &pac_encountered, pair<int, int> &lhs, set<pair<int, pair<int, int>>> &ec_params_and_vars, set<pair<int, pair<int, int>>> &ar_params_and_vars) const
 {
@@ -7564,16 +7573,11 @@ SecondDerivExternalFunctionNode::sameTefTermPredicate() const
 }
 
 VarExpectationNode::VarExpectationNode(DataTree &datatree_arg,
-                                       int symb_id_arg,
-                                       int forecast_horizon_arg,
-                                       const string &model_name_arg) :
+                                       string model_name_arg) :
   ExprNode(datatree_arg),
-  symb_id(symb_id_arg),
-  forecast_horizon(forecast_horizon_arg),
-  model_name(model_name_arg),
-  yidx(-1)
+  model_name{move(model_name_arg)}
 {
-  datatree.var_expectation_node_map[{ model_name, symb_id, forecast_horizon }] = this;
+  datatree.var_expectation_node_map[model_name] = this;
 }
 
 void
@@ -7581,7 +7585,8 @@ VarExpectationNode::computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReferenc
                                           map<NodeTreeReference, temporary_terms_t> &temp_terms_map,
                                           bool is_matlab, NodeTreeReference tr) const
 {
-  temp_terms_map[tr].insert(const_cast<VarExpectationNode *>(this));
+  cerr << "VarExpectationNode::computeTemporaryTerms not implemented." << endl;
+  exit(EXIT_FAILURE);
 }
 
 void
@@ -7592,22 +7597,21 @@ VarExpectationNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
                                           vector< vector<temporary_terms_t>> &v_temporary_terms,
                                           int equation) const
 {
-  expr_t this2 = const_cast<VarExpectationNode *>(this);
-  temporary_terms.insert(this2);
-  first_occurence[this2] = { Curr_block, equation };
-  v_temporary_terms[Curr_block][equation].insert(this2);
+  cerr << "VarExpectationNode::computeTemporaryTerms not implemented." << endl;
+  exit(EXIT_FAILURE);
 }
 
 expr_t
 VarExpectationNode::toStatic(DataTree &static_datatree) const
 {
-  return static_datatree.AddVariable(symb_id);
+  cerr << "VarExpectationNode::toStatic not implemented." << endl;
+  exit(EXIT_FAILURE);
 }
 
 expr_t
 VarExpectationNode::cloneDynamic(DataTree &dynamic_datatree) const
 {
-  return dynamic_datatree.AddVarExpectation(symb_id, forecast_horizon, model_name);
+  return dynamic_datatree.AddVarExpectation(model_name);
 }
 
 void
@@ -7620,101 +7624,110 @@ VarExpectationNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
 
   if (IS_LATEX(output_type))
     {
-      output << "VAR_" << model_name << LEFT_PAR(output_type)
-             << datatree.symbol_table.getTeXName(symb_id)
-             << "_{t+" << forecast_horizon << "}" << RIGHT_PAR(output_type);
+      output << "VAR_EXPECTATION(" << model_name << ')';
       return;
     }
 
-  if (checkIfTemporaryTermThenWrite(output, output_type, temporary_terms, temporary_terms_idxs))
-    return;
-
-  output << "dynamic_var_forecast_" << model_name << "_" << forecast_horizon << "(" << yidx + 1 << ")";
+  cerr << "VarExpectationNode::writeOutput not implemented for non-LaTeX." << endl;
+  exit(EXIT_FAILURE);
 }
 
 int
 VarExpectationNode::maxEndoLead() const
 {
-  return 0;
+  cerr << "VarExpectationNode::maxEndoLead not implemented." << endl;
+  exit(EXIT_FAILURE);
 }
 
 int
 VarExpectationNode::maxExoLead() const
 {
-  return 0;
+  cerr << "VarExpectationNode::maxExoLead not implemented." << endl;
+  exit(EXIT_FAILURE);
 }
 
 int
 VarExpectationNode::maxEndoLag() const
 {
-  return 0;
+  cerr << "VarExpectationNode::maxEndoLead not implemented." << endl;
+  exit(EXIT_FAILURE);
 }
 
 int
 VarExpectationNode::maxExoLag() const
 {
-  return 0;
+  cerr << "VarExpectationNode::maxExoLead not implemented." << endl;
+  exit(EXIT_FAILURE);
 }
 
 int
 VarExpectationNode::maxLead() const
 {
-  return 0;
+  cerr << "VarExpectationNode::maxLead not implemented." << endl;
+  exit(EXIT_FAILURE);
 }
 
 int
 VarExpectationNode::maxLag() const
 {
-  return 0;
+  cerr << "VarExpectationNode::maxLag not implemented." << endl;
+  exit(EXIT_FAILURE);
 }
 
 expr_t
 VarExpectationNode::undiff() const
 {
-  return const_cast<VarExpectationNode *>(this);
+  cerr << "VarExpectationNode::undiff not implemented." << endl;
+  exit(EXIT_FAILURE);
 }
 
 int
 VarExpectationNode::VarMinLag() const
 {
-  return 1;
+  cerr << "VarExpectationNode::VarMinLag not implemented." << endl;
+  exit(EXIT_FAILURE);
 }
 
 int
 VarExpectationNode::VarMaxLag(DataTree &static_datatree, set<expr_t> &static_lhs) const
 {
-  return 0;
+  cerr << "VarExpectationNode::VarMaxLag not implemented." << endl;
+  exit(EXIT_FAILURE);
 }
 
 int
 VarExpectationNode::PacMaxLag(vector<int> &lhs) const
 {
-  return 0;
+  cerr << "VarExpectationNode::PacMaxLag not implemented." << endl;
+  exit(EXIT_FAILURE);
 }
 
 expr_t
 VarExpectationNode::decreaseLeadsLags(int n) const
 {
-  return const_cast<VarExpectationNode *>(this);
+  cerr << "VarExpectationNode::decreaseLeadsLags not implemented." << endl;
+  exit(EXIT_FAILURE);
 }
 
 void
 VarExpectationNode::prepareForDerivation()
 {
-  preparedForDerivation = true;
-  // Come back
+  cerr << "VarExpectationNode::prepareForDerivation not implemented." << endl;
+  exit(EXIT_FAILURE);
 }
 
 expr_t
 VarExpectationNode::computeDerivative(int deriv_id)
 {
-  return datatree.Zero;
+  cerr << "VarExpectationNode::computeDerivative not implemented." << endl;
+  exit(EXIT_FAILURE);
 }
 
 expr_t
 VarExpectationNode::getChainRuleDerivative(int deriv_id, const map<int, expr_t> &recursive_variables)
 {
-  return datatree.Zero;
+  cerr << "VarExpectationNode::getChainRuleDerivative not implemented." << endl;
+  exit(EXIT_FAILURE);
 }
 
 bool
@@ -7726,17 +7739,14 @@ VarExpectationNode::containsExternalFunction() const
 double
 VarExpectationNode::eval(const eval_context_t &eval_context) const noexcept(false)
 {
-  auto it = eval_context.find(symb_id);
-  if (it == eval_context.end())
-    throw EvalException();
-
-  return it->second;
+  throw EvalException();
 }
 
 int
 VarExpectationNode::countDiffs() const
 {
-  return 0;
+  cerr << "VarExpectationNode::countDiffs not implemented." << endl;
+  exit(EXIT_FAILURE);
 }
 
 void
@@ -7759,9 +7769,8 @@ VarExpectationNode::collectDynamicVariables(SymbolType type_arg, set<pair<int, i
 void
 VarExpectationNode::collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const
 {
-  auto it = temporary_terms.find(const_cast<VarExpectationNode *>(this));
-  if (it != temporary_terms.end())
-    temporary_terms_inuse.insert(idx);
+  cerr << "VarExpectationNode::collectTemporary_terms not implemented." << endl;
+  exit(EXIT_FAILURE);
 }
 
 void
@@ -7777,31 +7786,36 @@ VarExpectationNode::compile(ostream &CompileCode, unsigned int &instruction_numb
 pair<int, expr_t >
 VarExpectationNode::normalizeEquation(int var_endo, vector<pair<int, pair<expr_t, expr_t>>> &List_of_Op_RHS) const
 {
-  return { 0, datatree.AddVariableInternal(symb_id, 0) };
+  cerr << "VarExpectationNode::normalizeEquation not implemented." << endl;
+  exit(EXIT_FAILURE);
 }
 
 expr_t
 VarExpectationNode::substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const
 {
-  return const_cast<VarExpectationNode *>(this);
+  cerr << "VarExpectationNode::substituteEndoLeadGreaterThanTwo not implemented." << endl;
+  exit(EXIT_FAILURE);
 }
 
 expr_t
 VarExpectationNode::substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
-  return const_cast<VarExpectationNode *>(this);
+  cerr << "VarExpectationNode::substituteEndoLagGreaterThanTwo not implemented." << endl;
+  exit(EXIT_FAILURE);
 }
 
 expr_t
 VarExpectationNode::substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const
 {
-  return const_cast<VarExpectationNode *>(this);
+  cerr << "VarExpectationNode::substituteExoLead not implemented." << endl;
+  exit(EXIT_FAILURE);
 }
 
 expr_t
 VarExpectationNode::substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
-  return const_cast<VarExpectationNode *>(this);
+  cerr << "VarExpectationNode::substituteExoLag not implemented." << endl;
+  exit(EXIT_FAILURE);
 }
 
 expr_t
@@ -7816,6 +7830,18 @@ VarExpectationNode::substituteAdl() const
   return const_cast<VarExpectationNode *>(this);
 }
 
+expr_t
+VarExpectationNode::substituteVarExpectation(const map<string, expr_t> &subst_table) const
+{
+  auto it = subst_table.find(model_name);
+  if (it == subst_table.end())
+    {
+      cerr << "ERROR: unknown model '" << model_name << "' used in var_expectation expression" << endl;
+      exit(EXIT_FAILURE);
+    }
+  return it->second;
+}
+
 void
 VarExpectationNode::findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const
 {
@@ -7848,7 +7874,8 @@ VarExpectationNode::substitutePacExpectation(map<const PacExpectationNode *, con
 expr_t
 VarExpectationNode::differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
-  return const_cast<VarExpectationNode *>(this);
+  cerr << "VarExpectationNode::differentiateForwardVars not implemented." << endl;
+  exit(EXIT_FAILURE);
 }
 
 bool
@@ -7860,13 +7887,15 @@ VarExpectationNode::containsPacExpectation() const
 bool
 VarExpectationNode::containsEndogenous() const
 {
-  return true;
+  cerr << "VarExpectationNode::containsEndogenous not implemented." << endl;
+  exit(EXIT_FAILURE);
 }
 
 bool
 VarExpectationNode::containsExogenous() const
 {
-  return false;
+  cerr << "VarExpectationNode::containsExogenous not implemented." << endl;
+  exit(EXIT_FAILURE);
 }
 
 bool
@@ -7878,7 +7907,8 @@ VarExpectationNode::isNumConstNodeEqualTo(double value) const
 expr_t
 VarExpectationNode::decreaseLeadsLagsPredeterminedVariables() const
 {
-  return const_cast<VarExpectationNode *>(this);
+  cerr << "VarExpectationNode::decreaseLeadsLagsPredeterminedVariables not implemented." << endl;
+  exit(EXIT_FAILURE);
 }
 
 bool
@@ -7890,31 +7920,37 @@ VarExpectationNode::isVariableNodeEqualTo(SymbolType type_arg, int variable_id,
 expr_t
 VarExpectationNode::replaceTrendVar() const
 {
-  return const_cast<VarExpectationNode *>(this);
+  cerr << "VarExpectationNode::replaceTrendVar not implemented." << endl;
+  exit(EXIT_FAILURE);
 }
 
 expr_t
 VarExpectationNode::detrend(int symb_id, bool log_trend, expr_t trend) const
 {
-  return const_cast<VarExpectationNode *>(this);
+  cerr << "VarExpectationNode::detrend not implemented." << endl;
+  exit(EXIT_FAILURE);
 }
 
 expr_t
 VarExpectationNode::removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const
 {
-  return const_cast<VarExpectationNode *>(this);
+  cerr << "VarExpectationNode::removeTrendLeadLag not implemented." << endl;
+  exit(EXIT_FAILURE);
 }
 
 bool
 VarExpectationNode::isInStaticForm() const
 {
-  return false;
+  cerr << "VarExpectationNode::isInStaticForm not implemented." << endl;
+  exit(EXIT_FAILURE);
 }
 
 bool
 VarExpectationNode::isVarModelReferenced(const string &model_info_name) const
 {
-  return model_name == model_info_name;
+  /* TODO: should check here whether the var_expectation_model is equal to the
+     argument; we probably need a VarModelTable class to do that elegantly */
+  return false;
 }
 
 void
@@ -7922,13 +7958,6 @@ VarExpectationNode::getEndosAndMaxLags(map<string, int> &model_endos_and_lags) c
 {
 }
 
-void
-VarExpectationNode::setVarExpectationIndex(map<string, pair<SymbolList, int>> &var_model_info)
-{
-  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();
-}
-
 void
 VarExpectationNode::walkPacParameters(bool &pac_encountered, pair<int, int> &lhs, set<pair<int, pair<int, int>>> &ec_params_and_vars, set<pair<int, pair<int, int>>> &ar_params_and_vars) const
 {
@@ -7956,12 +7985,7 @@ VarExpectationNode::writeJsonOutput(ostream &output,
                                     const deriv_node_temp_terms_t &tef_terms,
                                     const bool isdynamic) const
 {
-  output << "var_expectation("
-         << "forecast_horizon = " << forecast_horizon
-         << ", name = " << datatree.symbol_table.getName(symb_id)
-         << ", model_name = " << model_name
-         << ", yindex = " << yidx
-         << ")";
+  output << "var_expectation(" << model_name << ')';
 }
 
 PacExpectationNode::PacExpectationNode(DataTree &datatree_arg,
@@ -8271,6 +8295,12 @@ PacExpectationNode::substituteAdl() const
   return const_cast<PacExpectationNode *>(this);
 }
 
+expr_t
+PacExpectationNode::substituteVarExpectation(const map<string, expr_t> &subst_table) const
+{
+  return const_cast<PacExpectationNode *>(this);
+}
+
 void
 PacExpectationNode::findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const
 {
@@ -8371,11 +8401,6 @@ PacExpectationNode::getEndosAndMaxLags(map<string, int> &model_endos_and_lags) c
 {
 }
 
-void
-PacExpectationNode::setVarExpectationIndex(map<string, pair<SymbolList, int>> &var_model_info)
-{
-}
-
 expr_t
 PacExpectationNode::substituteStaticAuxiliaryVariable() const
 {
diff --git a/src/ExprNode.hh b/src/ExprNode.hh
index 6d595611ecd0d6982c1ba62f9b4388321495540a..e2710b7522ad39c007c39430d47c60aa789c9fdb 100644
--- a/src/ExprNode.hh
+++ b/src/ExprNode.hh
@@ -499,6 +499,9 @@ class ExprNode
       //! Substitute adl operator
       virtual expr_t substituteAdl() const = 0;
 
+      //! Substitute VarExpectation nodes
+      virtual expr_t substituteVarExpectation(const map<string, expr_t> &subst_table) const = 0;
+
       //! Substitute diff operator
       virtual void findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const = 0;
       virtual void findUnaryOpNodesForAuxVarCreation(DataTree &static_datatree, diff_table_t &nodes) const = 0;
@@ -520,9 +523,6 @@ class ExprNode
       //! 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, 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;
 
@@ -600,6 +600,7 @@ public:
   expr_t substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const override;
   expr_t substituteAdl() const override;
+  expr_t substituteVarExpectation(const map<string, expr_t> &subst_table) const override;
   void findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const override;
   void findUnaryOpNodesForAuxVarCreation(DataTree &static_datatree, diff_table_t &nodes) const override;
   expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
@@ -617,7 +618,6 @@ public:
   expr_t cloneDynamic(DataTree &dynamic_datatree) const override;
   expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const override;
   bool isInStaticForm() const override;
-  void setVarExpectationIndex(map<string, pair<SymbolList, int>> &var_model_info) override;
   void walkPacParameters(bool &pac_encountered, pair<int, int> &lhs, set<pair<int, pair<int, int>>> &ec_params_and_vars, set<pair<int, pair<int, int>>> &params_and_vars) const override;
   void addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg) override;
   void fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int equation_number_arg) override;
@@ -691,6 +691,7 @@ public:
   expr_t substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const override;
   expr_t substituteAdl() const override;
+  expr_t substituteVarExpectation(const map<string, expr_t> &subst_table) const override;
   void findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const override;
   void findUnaryOpNodesForAuxVarCreation(DataTree &static_datatree, diff_table_t &nodes) const override;
   expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
@@ -708,7 +709,6 @@ public:
   expr_t cloneDynamic(DataTree &dynamic_datatree) const override;
   expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const override;
   bool isInStaticForm() const override;
-  void setVarExpectationIndex(map<string, pair<SymbolList, int>> &var_model_info) override;
   void walkPacParameters(bool &pac_encountered, pair<int, int> &lhs, set<pair<int, pair<int, int>>> &ec_params_and_vars, set<pair<int, pair<int, int>>> &params_and_vars) const override;
   void addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg) override;
   void fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int equation_number_arg) override;
@@ -805,6 +805,7 @@ public:
   expr_t substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const override;
   expr_t substituteAdl() const override;
+  expr_t substituteVarExpectation(const map<string, expr_t> &subst_table) const override;
   void findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const override;
   bool createAuxVarForUnaryOpNode() const;
   void findUnaryOpNodesForAuxVarCreation(DataTree &static_datatree, diff_table_t &nodes) const override;
@@ -823,7 +824,6 @@ public:
   expr_t cloneDynamic(DataTree &dynamic_datatree) const override;
   expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const override;
   bool isInStaticForm() const override;
-  void setVarExpectationIndex(map<string, pair<SymbolList, int>> &var_model_info) override;
   void walkPacParameters(bool &pac_encountered, pair<int, int> &lhs, set<pair<int, pair<int, int>>> &ec_params_and_vars, set<pair<int, pair<int, int>>> &params_and_vars) const override;
   void addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg) override;
   void fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int equation_number_arg) override;
@@ -937,6 +937,7 @@ public:
   expr_t substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const override;
   expr_t substituteAdl() const override;
+  expr_t substituteVarExpectation(const map<string, expr_t> &subst_table) const override;
   void findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const override;
   void findUnaryOpNodesForAuxVarCreation(DataTree &static_datatree, diff_table_t &nodes) const override;
   expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
@@ -960,7 +961,6 @@ public:
   //! Returns the non-zero hand-side of an equation (that must have a hand side equal to zero)
   expr_t getNonZeroPartofEquation() const;
   bool isInStaticForm() const override;
-  void setVarExpectationIndex(map<string, pair<SymbolList, int>> &var_model_info) override;
   void walkPacParameters(bool &pac_encountered, pair<int, int> &lhs, set<pair<int, pair<int, int>>> &ec_params_and_vars, set<pair<int, pair<int, int>>> &params_and_vars) const override;
   void addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &ar_params_and_vars_arg) override;
   void fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int equation_number_arg) override;
@@ -1044,6 +1044,7 @@ public:
   expr_t substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const override;
   expr_t substituteAdl() const override;
+  expr_t substituteVarExpectation(const map<string, expr_t> &subst_table) const override;
   void findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const override;
   void findUnaryOpNodesForAuxVarCreation(DataTree &static_datatree, diff_table_t &nodes) const override;
   expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
@@ -1061,7 +1062,6 @@ public:
   expr_t cloneDynamic(DataTree &dynamic_datatree) const override;
   expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const override;
   bool isInStaticForm() const override;
-  void setVarExpectationIndex(map<string, pair<SymbolList, int>> &var_model_info) override;
   void walkPacParameters(bool &pac_encountered, pair<int, int> &lhs, set<pair<int, pair<int, int>>> &ec_params_and_vars, set<pair<int, pair<int, int>>> &params_and_vars) const override;
   void addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg) override;
   void fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int equation_number_arg) override;
@@ -1155,6 +1155,7 @@ public:
   expr_t substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const override;
   expr_t substituteAdl() const override;
+  expr_t substituteVarExpectation(const map<string, expr_t> &subst_table) const override;
   void findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const override;
   void findUnaryOpNodesForAuxVarCreation(DataTree &static_datatree, diff_table_t &nodes) const override;
   expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
@@ -1174,7 +1175,6 @@ public:
   expr_t cloneDynamic(DataTree &dynamic_datatree) const override = 0;
   expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const override;
   bool isInStaticForm() const override;
-  void setVarExpectationIndex(map<string, pair<SymbolList, int>> &var_model_info) override;
   void walkPacParameters(bool &pac_encountered, pair<int, int> &lhs, set<pair<int, pair<int, int>>> &ec_params_and_vars, set<pair<int, pair<int, int>>> &params_and_vars) const override;
   void addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg) override;
   void fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int equation_number_arg) override;
@@ -1312,12 +1312,9 @@ public:
 class VarExpectationNode : public ExprNode
 {
 private:
-  const int symb_id;
-  const int forecast_horizon;
-  const string &model_name;
-  int yidx;
+  const string model_name;
 public:
-  VarExpectationNode(DataTree &datatree_arg, int symb_id_arg, int forecast_horizon_arg, const string &model_name);
+  VarExpectationNode(DataTree &datatree_arg, string model_name_arg);
   void computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference>> &reference_count,
                                      map<NodeTreeReference, temporary_terms_t> &temp_terms_map,
                                      bool is_matlab, NodeTreeReference tr) const override;
@@ -1353,6 +1350,7 @@ public:
   expr_t substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const override;
   expr_t substituteAdl() const override;
+  expr_t substituteVarExpectation(const map<string, expr_t> &subst_table) const override;
   void findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const override;
   void findUnaryOpNodesForAuxVarCreation(DataTree &static_datatree, diff_table_t &nodes) const override;
   expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
@@ -1377,7 +1375,6 @@ public:
   expr_t detrend(int symb_id, bool log_trend, expr_t trend) const override;
   expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const override;
   bool isInStaticForm() const override;
-  void setVarExpectationIndex(map<string, pair<SymbolList, int>> &var_model_info) override;
   void walkPacParameters(bool &pac_encountered, pair<int, int> &lhs, set<pair<int, pair<int, int>>> &ec_params_and_vars, set<pair<int, pair<int, int>>> &params_and_vars) const override;
   void addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg) override;
   void fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int equation_number_arg) override;
@@ -1439,6 +1436,7 @@ public:
   expr_t substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const override;
   expr_t substituteAdl() const override;
+  expr_t substituteVarExpectation(const map<string, expr_t> &subst_table) const override;
   void findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const override;
   void findUnaryOpNodesForAuxVarCreation(DataTree &static_datatree, diff_table_t &nodes) const override;
   expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
@@ -1463,7 +1461,6 @@ public:
   expr_t detrend(int symb_id, bool log_trend, expr_t trend) const override;
   expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const override;
   bool isInStaticForm() const override;
-  void setVarExpectationIndex(map<string, pair<SymbolList, int>> &var_model_info) override;
   void walkPacParameters(bool &pac_encountered, pair<int, int> &lhs, set<pair<int, pair<int, int>>> &ec_params_and_vars, set<pair<int, pair<int, int>>> &params_and_vars) const override;
   void addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg) override;
   void fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int equation_number_arg) override;
diff --git a/src/ModFile.cc b/src/ModFile.cc
index b62db9d039f632a905addcbc379ef6abf5e22dd4..4f8d1e087906eaf8d8165799723f6293ecbca62b 100644
--- a/src/ModFile.cc
+++ b/src/ModFile.cc
@@ -444,11 +444,7 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const
      }
 
   if (!var_model_info_var_expectation.empty())
-    {
-      dynamic_model.setVarExpectationIndices(var_model_info_var_expectation);
-      dynamic_model.addEquationsForVar(var_model_info_var_expectation);
-    }
-  dynamic_model.fillVarExpectationFunctionsToWrite();
+    dynamic_model.addEquationsForVar(var_model_info_var_expectation);
 
   if (symbol_table.predeterminedNbr() > 0)
     dynamic_model.transformPredeterminedVariables();
@@ -511,6 +507,69 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const
         exit(EXIT_FAILURE);
     }
 
+  /* Handle var_expectation_model statements: collect information about them,
+     create the new corresponding parameters, and the expressions to replace
+     the var_expectation statements.
+     TODO: move information collection to checkPass(), within a new
+     VarModelTable class */
+  map<string, expr_t> var_expectation_subst_table;
+  for (auto & statement : statements)
+    {
+      auto vems = dynamic_cast<VarExpectationModelStatement *>(statement);
+      if (!vems)
+        continue;
+
+      auto &model_name = vems->model_name;
+
+      /* Find the corresponding VM statement and extract information for it;
+         ideally we should have a VarModelTable for that purpose */
+      VarModelStatement *vms = nullptr;
+      for (auto & statement2 : statements)
+        if ((vms = dynamic_cast<VarModelStatement *>(statement2))
+            && vms->name == vems->var_model_name)
+          break;
+      if (!vms)
+        {
+          cerr << "ERROR: var_expectation_model " << model_name << " refers to nonexistent " << vems->var_model_name << " var_model" << endl;
+          exit(EXIT_FAILURE);
+        }
+      /* The code below is duplicated further below in this function; but we
+         can't avoid that until the collecting of information about VARs is
+         moved *before* lead/lag substitutions (or even better, in checkPass()) */
+      vector<expr_t> lhs_expr_t;
+      vector<int> lhs, eqnumber;
+      vector<set<pair<int, int>>> rhs;
+      vector<bool> nonstationary;
+      vector<string> eqtags{var_model_eq_tags[var_model_name]};
+      dynamic_model.getVarModelVariablesFromEqTags(eqtags,
+                                                   eqnumber, lhs, lhs_expr_t, rhs, nonstationary);
+      int max_lag{original_model.getVarMaxLag(diff_static_model, eqnumber)};
+
+      /* Create auxiliary parameters and the expression to be substituted to
+         the var_expectations statement */
+      auto subst_expr = dynamic_model.Zero;
+      for (int lag = 1; lag <= max_lag; lag++)
+        for (auto variable : lhs)
+          {
+            string param_name = "var_expectation_model_" + model_name + '_' + symbol_table.getName(variable) + '_' + to_string(lag);
+            int new_param_id = symbol_table.addSymbol(param_name, SymbolType::parameter);
+            vems->aux_params_ids.push_back(new_param_id);
+
+            subst_expr = dynamic_model.AddPlus(subst_expr,
+                                               dynamic_model.AddTimes(dynamic_model.AddVariable(new_param_id),
+                                                                      dynamic_model.AddVariable(variable, -lag)));
+          }
+
+      if (var_expectation_subst_table.find(model_name) != var_expectation_subst_table.end())
+        {
+          cerr << "ERROR: model name '" << model_name << "' is used by several var_expectation_model statements" << endl;
+          exit(EXIT_FAILURE);
+        }
+      var_expectation_subst_table[model_name] = subst_expr;
+    }
+  // And finally perform the substitutions
+  dynamic_model.substituteVarExpectation(var_expectation_subst_table);
+
   if (mod_file_struct.stoch_simul_present
       || mod_file_struct.estimation_present
       || mod_file_struct.osr_present
@@ -532,6 +591,10 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const
       dynamic_model.substituteEndoLagGreaterThanTwo(true);
     }
 
+  /* TODO: most of this should rather be done in checkPass(). Also we should
+     probably add a VarModelTable class (similar to SymbolTable) for storing all
+     this information about VAR models, instead of putting it inside the
+     Statement(s) classes */
   for (auto & statement : statements)
     {
       auto *vms = dynamic_cast<VarModelStatement *>(statement);
@@ -997,11 +1060,6 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool clear_glo
       auto *lpass = dynamic_cast<LoadParamsAndSteadyStateStatement *>(statement);
       if (lpass && !no_static)
         static_model.writeAuxVarInitval(mOutputFile, oMatlabOutsideModel);
-
-      // Special treatement for Var Models
-      auto *vms = dynamic_cast<VarModelStatement *>(statement);
-      if (vms != nullptr)
-        vms->createVarModelMFunction(mOutputFile, dynamic_model.getVarExpectationFunctionsToWrite());
     }
 
   mOutputFile << "save('" << basename << "_results.mat', 'oo_', 'M_', 'options_');" << endl
diff --git a/src/ParsingDriver.cc b/src/ParsingDriver.cc
index 727a6c9d927ba56033c4a48789bdbebf3482a368..dc97c67c814596adee7641e9e6a3ef726ec8ea3f 100644
--- a/src/ParsingDriver.cc
+++ b/src/ParsingDriver.cc
@@ -2508,9 +2508,9 @@ ParsingDriver::add_expectation(const string &arg1, expr_t arg2)
 }
 
 expr_t
-ParsingDriver::add_var_expectation(const string &arg1, const string &arg2, const string &arg3)
+ParsingDriver::add_var_expectation(const string &model_name)
 {
-  return data_tree->AddVarExpectation(mod_file->symbol_table.getID(arg1), stoi(arg2), arg3);
+  return data_tree->AddVarExpectation(model_name);
 }
 
 expr_t
@@ -3228,3 +3228,45 @@ ParsingDriver::end_shock_groups(const string &name)
   mod_file->addStatement(new ShockGroupsStatement(shock_groups, name));
   shock_groups.clear();
 }
+
+void
+ParsingDriver::var_expectation_model()
+{
+  auto it = options_list.string_options.find("variable");
+  if (it == options_list.string_options.end())
+    error("You must pass the variable option to the var_expectation_model statement.");
+  auto variable = it->second;
+  check_symbol_is_endogenous(variable);
+
+  it = options_list.string_options.find("var_model_name");
+  if (it == options_list.string_options.end())
+    error("You must pass the var_model_name option to the var_expectation_model statement.");
+  auto var_model_name = it->second;
+
+  it = options_list.string_options.find("model_name");
+  if (it == options_list.string_options.end())
+    error("You must pass the model_name option to the var_expectation_model statement.");
+  auto model_name = it->second;
+
+  it = options_list.num_options.find("horizon");
+  if (it == options_list.num_options.end())
+    error("You must pass the horizon option to the var_expectation_model statement.");
+  auto horizon = it->second;
+
+  if (var_expectation_model_discount)
+    {
+      VariableNode *var;
+      if (!dynamic_cast<NumConstNode *>(var_expectation_model_discount)
+          && !((var = dynamic_cast<VariableNode *>(var_expectation_model_discount))
+               && var->get_type() == SymbolType::parameter))
+        error("The discount factor must be a constant expression or a parameter");
+    }
+  else
+    var_expectation_model_discount = data_tree->One;
+
+  mod_file->addStatement(new VarExpectationModelStatement(model_name, variable, var_model_name, horizon,
+                                                          var_expectation_model_discount, mod_file->symbol_table));
+
+  options_list.clear();
+  var_expectation_model_discount = nullptr;
+}
diff --git a/src/ParsingDriver.hh b/src/ParsingDriver.hh
index 808aad5cd59cfb9774c0cbd8114b150f57913702..e055199f5749d5fe44355c193f7dbccd8b04c244 100644
--- a/src/ParsingDriver.hh
+++ b/src/ParsingDriver.hh
@@ -297,6 +297,9 @@ public:
   map<pair<int, int>, double> covariance_number_restriction;
   map<pair<int, int>, pair<int, int>> covariance_pair_restriction;
 
+  //! Temporary storage for discount option of VAR_EXPECTATION_MODEL
+  expr_t var_expectation_model_discount{nullptr};
+
   //! Error handler with explicit location
   void error(const Dynare::parser::location_type &l, const string &m) __attribute__ ((noreturn));
   //! Error handler using saved location
@@ -688,8 +691,8 @@ public:
   expr_t add_power(expr_t arg1,  expr_t arg2);
   //! Writes token "E(arg1)(arg2)" to model tree
   expr_t add_expectation(const string &arg1,  expr_t arg2);
-  //! Writes token "VAR_EXPECTATION(arg1, arg2, arg3)" to model tree
-  expr_t add_var_expectation(const string &arg1, const string &arg2, const string &arg3);
+  //! Writes token "VAR_EXPECTATION(model_name)" to model tree
+  expr_t add_var_expectation(const string &model_name);
   //! Writes token "PAC_EXPECTATION(model_name, discount, growth)" to model tree
   expr_t add_pac_expectation(const string &var_model_name);
   //! Creates pac_model statement
@@ -836,6 +839,8 @@ public:
   void gmm_estimation();
   //! SMM Estimation statement
   void smm_estimation();
+  //! Add a var_expectation_model statement
+  void var_expectation_model();
 };
 
 #endif // ! PARSING_DRIVER_HH