diff --git a/src/ComputingTasks.cc b/src/ComputingTasks.cc
index 5947437de3ba17a5c7c628ea9b037c36a7a0b6df..cb94d7496655e8aeef83478ff525db36aaceaf67 100644
--- a/src/ComputingTasks.cc
+++ b/src/ComputingTasks.cc
@@ -4915,15 +4915,24 @@ GenerateIRFsStatement::writeJsonOutput(ostream &output) const
 }
 
 VarExpectationModelStatement::VarExpectationModelStatement(string model_name_arg,
-                                                           string variable_arg,
+                                                           expr_t expression_arg,
                                                            string aux_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)},
+  model_name{move(model_name_arg)}, expression{expression_arg},
   aux_model_name{move(aux_model_name_arg)}, horizon{move(horizon_arg)},
   discount{discount_arg}, symbol_table{symbol_table_arg}
 {
+  auto vpc = expression->matchLinearCombinationOfVariables();
+  for (const auto &it : vpc)
+    {
+      if (get<1>(it) != 0)
+        throw ExprNode::MatchFailureException{"lead/lags are not allowed"};
+      if (symbol_table.getType(get<0>(it)) != SymbolType::endogenous)
+        throw ExprNode::MatchFailureException{"Variable is not an endogenous"};
+      vars_params_constants.emplace_back(get<0>(it), get<2>(it), get<3>(it));
+    }
 }
 
 void
@@ -4931,9 +4940,28 @@ VarExpectationModelStatement::writeOutput(ostream &output, const string &basenam
 {
   string mstruct = "M_.var_expectation." + model_name;
   output << mstruct << ".auxiliary_model_name = '" << aux_model_name << "';" << endl
-         << mstruct << ".horizon = " << horizon << ';' << endl
-         << mstruct << ".variable = '" << variable << "';" << endl
-         << mstruct << ".variable_id = " << symbol_table.getTypeSpecificID(variable)+1 << ";" << endl;
+         << mstruct << ".horizon = " << horizon << ';' << endl;
+
+  ostringstream vars_list, params_list, constants_list;
+  for (auto it = vars_params_constants.begin(); it != vars_params_constants.end(); ++it)
+    {
+      if (it != vars_params_constants.begin())
+        {
+          vars_list << ", ";
+          params_list << ", ";
+          constants_list << ", ";
+        }
+      vars_list << symbol_table.getTypeSpecificID(get<0>(*it))+1;
+      if (get<1>(*it) == -1)
+        params_list << "NaN";
+      else
+        params_list << symbol_table.getTypeSpecificID(get<1>(*it))+1;
+      constants_list << get<2>(*it);
+    }
+  output << mstruct << ".expr.vars = [ " << vars_list.str() << " ];" << endl
+         << mstruct << ".expr.params = [ " << params_list.str() << " ];" << endl
+         << mstruct << ".expr.constants = [ " << constants_list.str() << " ];" << endl;
+
   auto disc_var = dynamic_cast<const VariableNode *>(discount);
   if (disc_var)
     output << mstruct << ".discount_index = " << symbol_table.getTypeSpecificID(disc_var->symb_id) + 1 << ';' << endl;
@@ -4954,7 +4982,9 @@ VarExpectationModelStatement::writeJsonOutput(ostream &output) const
 {
   output << "{\"statementName\": \"var_expectation_model\","
          << "\"model_name\": \"" << model_name << "\", "
-         << "\"variable\": \"" << variable << "\", "
+         << "\"expression\": \"";
+  expression->writeOutput(output);
+  output << "\", "
          << "\"auxiliary_model_name\": \"" << aux_model_name << "\", "
          << "\"horizon\": \"" << horizon << "\", "
          << "\"discount\": \"";
diff --git a/src/ComputingTasks.hh b/src/ComputingTasks.hh
index 09e52924231c8986e8a36ba32142a49c89b13f79..2afcbf766f8a0d606336058858e1804cf73952f9 100644
--- a/src/ComputingTasks.hh
+++ b/src/ComputingTasks.hh
@@ -1184,14 +1184,17 @@ public:
 class VarExpectationModelStatement : public Statement
 {
 public:
-  const string model_name, variable, aux_model_name, horizon;
+  const string model_name;
+  const expr_t expression;
+  const string aux_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
-
+private:
+  vector<tuple<int, int, double>> vars_params_constants;
 public:
-  VarExpectationModelStatement(string model_name_arg, string variable_arg, string aux_model_name_arg,
+  VarExpectationModelStatement(string model_name_arg, expr_t expression_arg, string aux_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;
diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 8cd9d89b6c4244a71aa78093300ab3822e891b67..ad9e6909511b80ad57c0b974de00ce880110f5d6 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -4271,7 +4271,16 @@ DynamicModel::walkPacParameters()
             {
               optim_share_index = *(optim_share.begin());
               optim_part->getPacOptimizingPart(lhs_orig_symb_id, ec_params_and_vars, ar_params_and_vars);
-              non_optim_vars_params_and_constants = non_optim_part->getPacNonOptimizingPart();
+              try
+                {
+                  non_optim_vars_params_and_constants = non_optim_part->matchLinearCombinationOfVariables();
+                }
+              catch (ExprNode::MatchFailureException &e)
+                {
+                  cerr << "Error in parsing non-optimizing agents part of PAC equation: "
+                       << e.message << endl;
+                  exit(EXIT_FAILURE);
+                }
             }
           equation->addParamInfoToPac(lhs,
                                       optim_share_index,
diff --git a/src/DynareBison.yy b/src/DynareBison.yy
index c42ea1178f418b9701361af434d7289bcefb5853..fac0a1b15f9eaefddf9ac0ed20a8bb1f2b0c1378 100644
--- a/src/DynareBison.yy
+++ b/src/DynareBison.yy
@@ -89,7 +89,7 @@ using tuple_4strings = tuple<string,string,string,string>;
 %token CALIB_SMOOTHER CHANGE_TYPE CHECK CONDITIONAL_FORECAST CONDITIONAL_FORECAST_PATHS CONF_SIG CONSTANT CONTROLLED_VAREXO CORR CUTOFF CYCLE_REDUCTION LOGARITHMIC_REDUCTION
 %token COMMA CONSIDER_ALL_ENDOGENOUS CONSIDER_ONLY_OBSERVED INITIAL_CONDITION_DECOMPOSITION
 %token DATAFILE FILE SERIES DET_COND_FORECAST DOUBLING DR_CYCLE_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_MAXITER DR_ALGO DROP DSAMPLE DYNASAVE DYNATYPE CALIBRATION DIFFERENTIATE_FORWARD_VARS
-%token END ENDVAL EQUAL ESTIMATION ESTIMATED_PARAMS ESTIMATED_PARAMS_BOUNDS ESTIMATED_PARAMS_INIT EXTENDED_PATH ENDOGENOUS_PRIOR
+%token END ENDVAL EQUAL ESTIMATION ESTIMATED_PARAMS ESTIMATED_PARAMS_BOUNDS ESTIMATED_PARAMS_INIT EXTENDED_PATH ENDOGENOUS_PRIOR EXPRESSION
 %token FILENAME DIRNAME FILTER_STEP_AHEAD FILTERED_VARS FIRST_OBS LAST_OBS SET_TIME OSR_PARAMS_BOUNDS KEEP_KALMAN_ALGO_IF_SINGULARITY_IS_DETECTED
 %token <string> FLOAT_NUMBER DATES
 %token DEFAULT FIXED_POINT OPT_ALGO
@@ -421,6 +421,8 @@ var_expectation_model_options_list : var_expectation_model_option
 
 var_expectation_model_option : VARIABLE EQUAL symbol
                                { driver.option_str("variable", $3); }
+                             | EXPRESSION EQUAL expression
+                               { driver.var_expectation_model_expression = $3; }
                              | AUXILIARY_MODEL_NAME EQUAL symbol
                                { driver.option_str("auxiliary_model_name", $3); }
                              | HORIZON EQUAL INT_NUMBER
diff --git a/src/DynareFlex.ll b/src/DynareFlex.ll
index 7c2bce41e62ea408dd3b5a9f14a6a202fcb98eac..186889af93e96be9195804bcbea0412e3b9ba0a0 100644
--- a/src/DynareFlex.ll
+++ b/src/DynareFlex.ll
@@ -431,6 +431,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>no_posterior_kernel_density {return token::NO_POSTERIOR_KERNEL_DENSITY;}
 <DYNARE_STATEMENT>rescale_prediction_error_covariance {return token::RESCALE_PREDICTION_ERROR_COVARIANCE;}
 <DYNARE_STATEMENT>use_penalized_objective_for_hessian {return token::USE_PENALIZED_OBJECTIVE_FOR_HESSIAN;}
+<DYNARE_STATEMENT>expression    {return token::EXPRESSION;}
 
 <DYNARE_STATEMENT>alpha {
   yylval->build<string>(yytext);
diff --git a/src/ExprNode.cc b/src/ExprNode.cc
index f1ddd2c54e79f8baba96472f55f57134e8f8c389..a415b7f9e7e57615d34985d2c1705cb7d1363d7f 100644
--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
@@ -9640,7 +9640,7 @@ BinaryOpNode::matchVTCTPHelper(int &var_id, int &lag, int &param_id, double &con
 }
 
 vector<tuple<int, int, int, double>>
-ExprNode::getPacNonOptimizingPart() const
+ExprNode::matchLinearCombinationOfVariables() const
 {
   vector<pair<expr_t, int>> terms;
   decomposeAdditiveTerms(terms);
@@ -9648,19 +9648,12 @@ ExprNode::getPacNonOptimizingPart() const
   vector<tuple<int, int, int, double>> result;
 
   for (const auto &it : terms)
-    try
-      {
-        expr_t term = it.first;
-        int sign = it.second;
-        auto m = term->matchVariableTimesConstantTimesParam();
-        get<3>(m) *= sign;
-        result.push_back(m);
-      }
-    catch (MatchFailureException &e)
-      {
-        cerr << "ExprNode::getPacNonOptimizingPart: Error in parsing PAC equation: "
-             << e.message << endl;
-        exit(EXIT_FAILURE);
-      }
-    return result;
+    {
+      expr_t term = it.first;
+      int sign = it.second;
+      auto m = term->matchVariableTimesConstantTimesParam();
+      get<3>(m) *= sign;
+      result.push_back(m);
+    }
+  return result;
 }
diff --git a/src/ExprNode.hh b/src/ExprNode.hh
index 5c541bf2148ad8ff07bda875336c9c99c143a9e6..3927c6d4cc63a1e960e31d43bb58079e79b24896 100644
--- a/src/ExprNode.hh
+++ b/src/ExprNode.hh
@@ -560,11 +560,12 @@ class ExprNode
       virtual void getPacOptimizingPart(int lhs_orig_symb_id, pair<int, pair<vector<int>, vector<bool>>> &ec_params_and_vars,
                                         set<pair<int, pair<int, int>>> &params_and_vars) const = 0;
 
-      //! Analyzes the non optimizing part of PAC equation
+      //! Matches a linear combination of variables, where scalars can be constant*parameter
       /*! Returns a list of (variable_id, lag, param_id, constant)
           corresponding to the terms in the expression. When there is no
-          parameter in a term, param_id == -1 */
-      vector<tuple<int, int, int, double>> getPacNonOptimizingPart() const;
+          parameter in a term, param_id == -1.
+          Can throw a MatchFailureException. */
+      vector<tuple<int, int, int, double>> matchLinearCombinationOfVariables() const;
       //! Returns true if expression is of the form:
       //! param * (endog op endog op ...) + param * (endog op endog op ...) + ...
       virtual bool isParamTimesEndogExpr() const = 0;
diff --git a/src/ParsingDriver.cc b/src/ParsingDriver.cc
index c5c4c56fbbc126ac338fb1ef76115feb237c202c..482e5f90bd31cc12c2f02bc819e4c3f54516faa9 100644
--- a/src/ParsingDriver.cc
+++ b/src/ParsingDriver.cc
@@ -3348,10 +3348,14 @@ 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);
+  if (it == options_list.string_options.end() && !var_expectation_model_expression)
+    error("You must pass either the 'variable' or the 'expression' option to the var_expectation_model statement.");
+  if (it != options_list.string_options.end())
+    {
+      if (var_expectation_model_expression)
+        error("You can't pass both the 'variable' or the 'expression' options to the var_expectation_model statement.");
+      var_expectation_model_expression = data_tree->AddVariable(mod_file->symbol_table.getID(it->second));
+    }
 
   it = options_list.string_options.find("auxiliary_model_name");
   if (it == options_list.string_options.end())
@@ -3379,9 +3383,18 @@ ParsingDriver::var_expectation_model()
   else
     var_expectation_model_discount = data_tree->One;
 
-  mod_file->addStatement(make_unique<VarExpectationModelStatement>(model_name, variable, var_model_name, horizon,
-                                                                   var_expectation_model_discount, mod_file->symbol_table));
+  try
+    {
+      mod_file->addStatement(make_unique<VarExpectationModelStatement>(model_name, var_expectation_model_expression,
+                                                                       var_model_name, horizon,
+                                                                       var_expectation_model_discount, mod_file->symbol_table));
+    }
+  catch (ExprNode::MatchFailureException &e)
+    {
+      error("expression in var_expectation_model is not of the expected form: " + e.message);
+    }
 
   options_list.clear();
   var_expectation_model_discount = nullptr;
+  var_expectation_model_expression = nullptr;
 }
diff --git a/src/ParsingDriver.hh b/src/ParsingDriver.hh
index c8fb1969602fd19ac29b21e8efc5bc64bf97b857..39648a7a9a5220c9a4e65c10dc365287a8b20c8d 100644
--- a/src/ParsingDriver.hh
+++ b/src/ParsingDriver.hh
@@ -309,6 +309,8 @@ public:
   map<pair<int, int>, double> covariance_number_restriction;
   map<pair<int, int>, pair<int, int>> covariance_pair_restriction;
 
+  //! Temporary storage for "expression" option of VAR_EXPECTATION_MODEL
+  expr_t var_expectation_model_expression{nullptr};
   //! Temporary storage for discount option of VAR_EXPECTATION_MODEL
   expr_t var_expectation_model_discount{nullptr};