diff --git a/src/CodeInterpreter.hh b/src/CodeInterpreter.hh
index ad9149797e4efea110d8e5ce15b82cfa7b5e0a8e..d7f533b80b62639209bfb4de8715c85d7c73de18 100644
--- a/src/CodeInterpreter.hh
+++ b/src/CodeInterpreter.hh
@@ -137,7 +137,7 @@ enum class SymbolType
    // Value 17 is unused for the time being (but could be reused)
 
    epilogue = 18, //!< Variables created in epilogue block
-   excludedVariable = 19 //!< Type to use when an equation is excluded via include/exclude_eqs and the LHS variable is not used elsewhere in the model
+   excludedVariable = 19 //!< Variable excluded via model_remove/include_eqs/exclude_eqs
   };
 
 enum class ExpressionType
diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 0f76199caeb4e289d0b6d82251f1fde62ef104ee..39a6ce520ac9c57c3b25e242f8b1791364b1853a 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -2356,18 +2356,23 @@ DynamicModel::writeDynamicJacobianNonZeroElts(const string &basename) const
   output.close();
 }
 
-void
-DynamicModel::parseIncludeExcludeEquations(const string &inc_exc_eq_tags,
-                                           set<pair<string, string>> &eq_tag_set, bool exclude_eqs)
+vector<pair<string, string>>
+DynamicModel::parseIncludeExcludeEquations(const string &inc_exc_option_value, bool exclude_eqs)
 {
+  auto removeLeadingTrailingWhitespace = [](string &str)
+  {
+    str.erase(0, str.find_first_not_of("\t\n\v\f\r "));
+    str.erase(str.find_last_not_of("\t\n\v\f\r ") + 1);
+  };
+
   string tags;
-  if (filesystem::exists(inc_exc_eq_tags))
+  if (filesystem::exists(inc_exc_option_value))
     {
       ifstream exclude_file;
-      exclude_file.open(inc_exc_eq_tags, ifstream::in);
+      exclude_file.open(inc_exc_option_value, ifstream::in);
       if (!exclude_file.is_open())
         {
-          cerr << "ERROR: Could not open " << inc_exc_eq_tags << endl;
+          cerr << "ERROR: Could not open " << inc_exc_option_value << endl;
           exit(EXIT_FAILURE);
         }
 
@@ -2397,12 +2402,12 @@ DynamicModel::parseIncludeExcludeEquations(const string &inc_exc_eq_tags,
         }
     }
   else
-    tags = inc_exc_eq_tags;
+    tags = inc_exc_option_value;
   removeLeadingTrailingWhitespace(tags);
 
   if (tags.front() == '[' && tags.back() != ']')
     {
-      cerr << "Error: " << (exclude_eqs ? "exclude_eqs" : "include_eqs")
+      cerr << "ERROR: " << (exclude_eqs ? "exclude_eqs" : "include_eqs")
            << ": if the first character is '[' the last must be ']'" << endl;
       exit(EXIT_FAILURE);
     }
@@ -2433,11 +2438,12 @@ DynamicModel::parseIncludeExcludeEquations(const string &inc_exc_eq_tags,
   regex r(R"((\s*)" + quote_regex + "|" + non_quote_regex + R"(\s*)(,\s*()" + quote_regex + "|" + non_quote_regex + R"()\s*)*)");
   if (!regex_match(tags, r))
     {
-      cerr << "Error: " << (exclude_eqs ? "exclude_eqs" : "include_eqs")
+      cerr << "ERROR: " << (exclude_eqs ? "exclude_eqs" : "include_eqs")
            << ": argument is of incorrect format." << endl;
       exit(EXIT_FAILURE);
     }
 
+  vector<pair<string, string>> eq_tag_set;
   regex s(quote_regex + "|" + non_quote_regex);
   for (auto it = sregex_iterator(tags.begin(), tags.end(), s);
        it != sregex_iterator(); ++it)
@@ -2445,66 +2451,171 @@ DynamicModel::parseIncludeExcludeEquations(const string &inc_exc_eq_tags,
       auto str = it->str();
       if (str[0] == '\'' && str[str.size()-1] == '\'')
         str = str.substr(1, str.size()-2);
-      eq_tag_set.insert({tagname, str});
+      eq_tag_set.emplace_back(tagname, str);
     }
+  return eq_tag_set;
 }
 
-void
-DynamicModel::includeExcludeEquations(const string &eqs, bool exclude_eqs)
+vector<int>
+DynamicModel::removeEquationsHelper(set<pair<string, string>> &listed_eqs_by_tag, bool exclude_eqs,
+                                    bool excluded_vars_change_type,
+                                    vector<BinaryOpNode *> &all_equations,
+                                    vector<int> &all_equations_lineno,
+                                    EquationTags &all_equation_tags, bool static_equations) const
 {
-  if (eqs.empty())
-    return;
+  if (all_equations.empty())
+    return {};
+
+  /* Try to convert the list of equations by tags into a list of equation
+     numbers.
+     The tag pairs that match an equation are removed from the list, so that
+     the caller knows which tag pairs have not been handled. */
+  set<int> listed_eqs_by_number;
+  for (auto it = listed_eqs_by_tag.begin(); it != listed_eqs_by_tag.end();)
+    if (auto tmp = all_equation_tags.getEqnsByTag(it->first, it->second);
+        !tmp.empty())
+      {
+        listed_eqs_by_number.insert(tmp.begin(), tmp.end());
+        it = listed_eqs_by_tag.erase(it);
+      }
+    else
+      ++it;
+
+  // Compute the indices of equations to be actually deleted
+  set<int> eqs_to_delete_by_number;
+  if (exclude_eqs)
+    eqs_to_delete_by_number = listed_eqs_by_number;
+  else
+    for (size_t i = 0; i < all_equations.size(); i++)
+      if (listed_eqs_by_number.find(i) == listed_eqs_by_number.end())
+        eqs_to_delete_by_number.insert(i);
+
+  // remove from equations, equations_lineno, equation_tags
+  vector<BinaryOpNode *> new_equations;
+  vector<int> new_equations_lineno;
+  map<int, int> old_eqn_num_2_new;
+  vector<int> excluded_vars;
+  for (size_t i = 0; i < all_equations.size(); i++)
+    if (eqs_to_delete_by_number.find(i) != eqs_to_delete_by_number.end())
+      {
+        if (excluded_vars_change_type)
+          {
+            if (auto tmp = all_equation_tags.getTagValueByEqnAndKey(i, "endogenous"); !tmp.empty())
+              excluded_vars.push_back(symbol_table.getID(tmp));
+            else
+              {
+                set<int> result;
+                all_equations[i]->arg1->collectVariables(SymbolType::endogenous, result);
+                if (result.size() == 1)
+                  excluded_vars.push_back(*result.begin());
+                else
+                  {
+                    cerr << "ERROR: Equation " << i+1
+                         << " has been excluded but it does not have a single variable on its left-hand side or an `endogenous` tag" << endl;
+                    exit(EXIT_FAILURE);
+                  }
+              }
+          }
+      }
+    else
+      {
+        new_equations.emplace_back(all_equations[i]);
+        old_eqn_num_2_new[i] = new_equations.size() - 1;
+        new_equations_lineno.emplace_back(all_equations_lineno[i]);
+      }
+  int n_excl = all_equations.size() - new_equations.size();
+
+  all_equations = new_equations;
+  all_equations_lineno = new_equations_lineno;
+
+  all_equation_tags.erase(eqs_to_delete_by_number, old_eqn_num_2_new);
+
+  if (!static_equations)
+    for (size_t i = 0; i < excluded_vars.size(); i++)
+      for (size_t j = i+1; j < excluded_vars.size(); j++)
+        if (excluded_vars[i] == excluded_vars[j])
+          {
+            cerr << "ERROR: Variable " << symbol_table.getName(i) << " was excluded twice"
+                 << " via a model_remove or model_replace statement, or via the include_eqs or exclude_eqs option" << endl;
+            exit(EXIT_FAILURE);
+          }
+
+  cout << "Excluded " << n_excl << (static_equations ? " static " : " dynamic ")
+       << "equation" << (n_excl > 1 ? "s" : "") << " via model_remove or model_replace statement, or via include_eqs or exclude_eqs option" << endl;
 
-  set<pair<string, string>> eq_tag_set;
-  parseIncludeExcludeEquations(eqs, eq_tag_set, exclude_eqs);
+  return excluded_vars;
+}
+
+void
+DynamicModel::removeEquations(const vector<pair<string, string>> &listed_eqs_by_tag, bool exclude_eqs,
+                              bool excluded_vars_change_type)
+{
+  /* Convert the const vector to a (mutable) set */
+  set<pair<string, string>> listed_eqs_by_tag2;
+  copy(listed_eqs_by_tag.begin(), listed_eqs_by_tag.end(), inserter(listed_eqs_by_tag2, listed_eqs_by_tag2.end()));
 
-  vector<int> excluded_vars
-    = ModelTree::includeExcludeEquations(eq_tag_set, exclude_eqs,
-                                         equations, equations_lineno,
-                                         equation_tags, false);
+  vector<int> excluded_vars = removeEquationsHelper(listed_eqs_by_tag2, exclude_eqs,
+                                                    excluded_vars_change_type,
+                                                    equations, equations_lineno,
+                                                    equation_tags, false);
 
   // Ignore output because variables are not excluded when equations marked 'static' are excluded
-  ModelTree::includeExcludeEquations(eq_tag_set, exclude_eqs,
-                                     static_only_equations, static_only_equations_lineno,
-                                     static_only_equations_equation_tags, true);
+  removeEquationsHelper(listed_eqs_by_tag2, exclude_eqs, excluded_vars_change_type,
+                        static_only_equations, static_only_equations_lineno,
+                        static_only_equations_equation_tags, true);
 
-  if (!eq_tag_set.empty())
+  if (!listed_eqs_by_tag2.empty())
     {
-      cerr << "ERROR: " << (exclude_eqs ? "exclude_eqs" : "include_eqs") << ": The equations specified by `";
-      cerr << eq_tag_set.begin()->first << "= ";
-      for (auto &it : eq_tag_set)
-        cerr << it.second << ", ";
-      cerr << "` were not found." << endl;
+      cerr << "ERROR: model_remove/model_replace/exclude_eqs/include_eqs: The equations specified by" << endl;
+      for (const auto &[tagname, tagvalue] : listed_eqs_by_tag)
+        cerr << " " << tagname << "=" << tagvalue << endl;
+      cerr << "were not found." << endl;
       exit(EXIT_FAILURE);
     }
 
-  if (staticOnlyEquationsNbr() != dynamicOnlyEquationsNbr())
+  if (excluded_vars_change_type)
     {
-      cerr << "ERROR: " << (exclude_eqs ? "exclude_eqs" : "include_eqs")
-           << ": You must remove the same number of equations marked `static` as equations marked `dynamic`." << endl;
-      exit(EXIT_FAILURE);
-    }
+      // Collect list of used variables in updated list of equations
+      set<int> eqn_vars;
+      for (auto eqn : equations)
+        eqn->collectVariables(SymbolType::endogenous, eqn_vars);
+      for (auto eqn : static_only_equations)
+        eqn->collectVariables(SymbolType::endogenous, eqn_vars);
 
-  // Collect list of used variables in updated list of equations
-  set<pair<int, int>> eqn_vars;
-  for (const auto &eqn : equations)
-    eqn->collectDynamicVariables(SymbolType::endogenous, eqn_vars);
-  for (const auto &eqn : static_only_equations)
-    eqn->collectDynamicVariables(SymbolType::endogenous, eqn_vars);
-
-  // Change LHS variable type of excluded equation if it is used in an eqution that has been kept
-  for (auto ev : excluded_vars)
-    {
-      bool found = false;
-      for (const auto &it : eqn_vars)
-        if (it.first == ev)
+      /* Change type of endogenous variables determined by excluded equations.
+         They become exogenous if they are still used somewhere, otherwise they are
+         completely excluded from the model. */
+      for (auto ev : excluded_vars)
+        if (eqn_vars.find(ev) != eqn_vars.end())
           {
             symbol_table.changeType(ev, SymbolType::exogenous);
-            found = true;
-            break;
+            cerr << "Variable '" << symbol_table.getName(ev) << "' turned into an exogenous, as its defining equation has been removed (but it still appears in an equation)" << endl;
           }
-      if (!found)
-        symbol_table.changeType(ev, SymbolType::excludedVariable);
+        else
+          {
+            symbol_table.changeType(ev, SymbolType::excludedVariable);
+            cerr << "Variable '" << symbol_table.getName(ev) << "' has been excluded from the model, as its defining equation has been removed and it appears nowhere else" << endl;
+          }
+    }
+}
+
+void
+DynamicModel::includeExcludeEquations(const string &inc_exc_option_value, bool exclude_eqs)
+{
+  if (inc_exc_option_value.empty())
+    return;
+
+  auto listed_eqs_by_tag = parseIncludeExcludeEquations(inc_exc_option_value, exclude_eqs);
+
+  removeEquations(listed_eqs_by_tag, exclude_eqs, true);
+
+  /* There is already a check about #static and #dynamic in
+     ModFile::checkPass(), but the present method is called from
+     ModFile::transformPass(), so we must do the check again */
+  if (staticOnlyEquationsNbr() != dynamicOnlyEquationsNbr())
+    {
+      cerr << "ERROR: exclude_eqs/include_eqs: You must remove the same number of equations marked `static` as equations marked `dynamic`." << endl;
+      exit(EXIT_FAILURE);
     }
 }
 
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index 09d4d90086d06f2aa3955715b6b274f7f6a9a64d..05584e4270bb6d81ffb5480277606df1f5f72a43 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -244,32 +244,52 @@ private:
       pointers into their equivalent in the new tree */
   void copyHelper(const DynamicModel &m);
 
-  // Internal helper functions for includeExcludeEquations()
-  /*! Handles parsing of argument passed to exclude_eqs/include_eqs*/
-  /*
-    Expects command line arguments of the form:
+  /* Handles parsing of argument passed to exclude_eqs/include_eqs.
+
+    The argument inc_exc_option_value should be of one of the following forms:
       * filename.txt
       * eq1
       * ['eq 1', 'eq 2']
       * [tagname='eq 1']
       * [tagname=('eq 1', 'eq 2')]
-    If argument is a file, the file should be formatted as:
+    If argument is a filename, the file should be formatted as:
         eq 1
         eq 2
     OR
         tagname=
         X
         Y
-   */
-  void parseIncludeExcludeEquations(const string &inc_exc_eq_tags, set<pair<string, string>> &eq_tag_set, bool exclude_eqs);
 
-  // General function that removes leading/trailing whitespace from a string
-  inline void
-  removeLeadingTrailingWhitespace(string &str)
-  {
-    str.erase(0, str.find_first_not_of("\t\n\v\f\r "));
-    str.erase(str.find_last_not_of("\t\n\v\f\r ") + 1);
-  }
+    The boolean exclude_eqs should be true if we are in the exclude_eqs case,
+    false in the include_eqs case (this only affects error messages).
+
+    Returns a set of pairs (tag name, tag value) corresponding to the set of
+    equations to be included or excluded.
+   */
+  static vector<pair<string, string>> parseIncludeExcludeEquations(const string &inc_exc_option_value, bool exclude_eqs);
+
+  /* Helper for the removeEquations() method.
+     listed_eqs_by_tag is the list of (tag name, tag value) pairs corresponding
+     to the option value, exclude_eqs is a boolean indicating whether we’re
+     excluding or including, and excluded_vars_change_type is a boolean
+     indicating whether to compute variables to be excluded.
+
+     The all_equations* arguments will be modified by the routine by excluding
+     equations. They are either the main structures for storing equations in
+     ModelTree, or their counterpart for static-only equations. The
+     static_equations boolean indicates when we are in the latter case.
+     The listed_eqs_by_tag structure will be updated by removing those tag
+     pairs that have been matched with equations in the all_equations*
+     argument*.
+
+     Returns a list of excluded variables (empty if
+     excluded_vars_change_type=false) */
+  vector<int> removeEquationsHelper(set<pair<string, string>> &listed_eqs_by_tag,
+                                    bool exclude_eqs, bool excluded_vars_change_type,
+                                    vector<BinaryOpNode *> &all_equations,
+                                    vector<int> &all_equations_lineno,
+                                    EquationTags &all_equation_tags,
+                                    bool static_equations) const;
 
   //! Compute autoregressive matrices of trend component models
   /* The algorithm uses matching rules over expression trees. It cannot handle
@@ -409,8 +429,14 @@ public:
   //! Set the max leads/lags of the original model
   void setLeadsLagsOrig();
 
-  //! Removes equations from the model according to name tags
-  void includeExcludeEquations(const string &eqs, bool exclude_eqs);
+  //! Implements the include_eqs/exclude_eqs options
+  void includeExcludeEquations(const string &inc_exc_option_value, bool exclude_eqs);
+
+  /* Removes equations from the model (identified by their name tags).
+     Used for include_eqs/exclude_eqs options and for model_remove and
+     model_replace blocks */
+  void removeEquations(const vector<pair<string, string>> &listed_eqs_by_tag, bool exclude_eqs,
+                       bool excluded_vars_change_type);
 
   //! Replaces model equations with derivatives of Lagrangian w.r.t. endogenous
   void computeRamseyPolicyFOCs(const StaticModel &static_model);
diff --git a/src/DynareBison.yy b/src/DynareBison.yy
index 8e42c8aa8579864bc3a4b0489c1919f2ba4d4817..dc1db9dc379fe168129a2f0f91a1c2c06b948f10 100644
--- a/src/DynareBison.yy
+++ b/src/DynareBison.yy
@@ -183,7 +183,7 @@ class ParsingDriver;
 %token NO_IDENTIFICATION_STRENGTH NO_IDENTIFICATION_REDUCEDFORM NO_IDENTIFICATION_MOMENTS
 %token NO_IDENTIFICATION_MINIMAL NO_IDENTIFICATION_SPECTRUM NORMALIZE_JACOBIANS GRID_NBR
 %token TOL_RANK TOL_DERIV TOL_SV CHECKS_VIA_SUBSETS MAX_DIM_SUBSETS_GROUPS ZERO_MOMENTS_TOLERANCE
-%token MAX_NROWS SQUEEZE_SHOCK_DECOMPOSITION WITH_EPILOGUE
+%token MAX_NROWS SQUEEZE_SHOCK_DECOMPOSITION WITH_EPILOGUE MODEL_REMOVE MODEL_REPLACE
 
 %token <vector<string>> SYMBOL_VEC
 
@@ -205,7 +205,7 @@ class ParsingDriver;
 %type <PriorDistributions> prior_pdf prior_distribution
 %type <pair<expr_t,expr_t>> calibration_range
 %type <pair<string,string>> named_var_elem subsamples_eq_opt integer_range_w_inf
-%type <vector<pair<string,string>>> named_var named_var_1
+%type <vector<pair<string,string>>> named_var named_var_1 tag_pair_list_for_selection
 %type <tuple<string,string,string,string>> prior_eq_opt options_eq_opt
 %type <vector<pair<int, int>>> period_list
 %type <vector<expr_t>> matched_moments_list value_list
@@ -342,6 +342,8 @@ statement : parameters
           | compilation_setup
           | matched_moments
           | occbin_constraints
+          | model_remove
+          | model_replace
           ;
 
 dsample : DSAMPLE INT_NUMBER ';'
@@ -1109,6 +1111,30 @@ comma_hand_side : hand_side
 pound_expression: '#' symbol EQUAL hand_side ';'
                   { driver.declare_and_init_model_local_variable($2, $4); };
 
+model_remove : MODEL_REMOVE '(' tag_pair_list_for_selection ')' ';'
+               { driver.model_remove($3); };
+
+model_replace : MODEL_REPLACE '(' tag_pair_list_for_selection ')' ';'
+               { driver.begin_model_replace($3); }
+               equation_list END ';'
+               { driver.end_model(); };
+
+tag_pair_list_for_selection : QUOTED_STRING
+                              { $$ = { { "name", $1 } }; }
+                            | symbol EQUAL QUOTED_STRING
+                              { $$ = { { $1, $3 } }; }
+                            | tag_pair_list_for_selection COMMA QUOTED_STRING
+                              {
+                                $$ = $1;
+                                $$.emplace_back("name", $3);
+                              }
+                            | tag_pair_list_for_selection COMMA symbol EQUAL QUOTED_STRING
+                              {
+                                $$ = $1;
+                                $$.emplace_back($3, $5);
+                              }
+                            ;
+
 shocks : SHOCKS ';' shock_list END ';' { driver.end_shocks(false); }
        | SHOCKS '(' OVERWRITE ')' ';' shock_list END ';' { driver.end_shocks(true); }
        | SHOCKS '(' OVERWRITE ')' ';'  END ';' { driver.end_shocks(true); }
diff --git a/src/DynareFlex.ll b/src/DynareFlex.ll
index e19aa7ea8d80c4fda784b7b6b50704278ca21f18..3150bfa6aafd3d254dd937c81452bd43f15c9911 100644
--- a/src/DynareFlex.ll
+++ b/src/DynareFlex.ll
@@ -192,6 +192,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <INITIAL>perfect_foresight_with_expectation_errors_setup {BEGIN DYNARE_STATEMENT; return token::PERFECT_FORESIGHT_WITH_EXPECTATION_ERRORS_SETUP;}
 <INITIAL>perfect_foresight_with_expectation_errors_solver {BEGIN DYNARE_STATEMENT; return token::PERFECT_FORESIGHT_WITH_EXPECTATION_ERRORS_SOLVER;}
 <INITIAL>compilation_setup {BEGIN DYNARE_STATEMENT; return token::COMPILATION_SETUP;}
+<INITIAL>model_remove {BEGIN DYNARE_STATEMENT; return token::MODEL_REMOVE;}
 
 <DYNARE_STATEMENT>; {
   if (!sigma_e)
@@ -231,6 +232,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <INITIAL>generate_irfs {BEGIN DYNARE_BLOCK; return token::GENERATE_IRFS;}
 <INITIAL>matched_moments {BEGIN DYNARE_BLOCK; return token::MATCHED_MOMENTS;}
 <INITIAL>occbin_constraints {BEGIN DYNARE_BLOCK; return token::OCCBIN_CONSTRAINTS;}
+<INITIAL>model_replace {BEGIN DYNARE_BLOCK; return token::MODEL_REPLACE;}
 
  /* For the semicolon after an "end" keyword */
 <INITIAL>; {return Dynare::parser::token_type (yytext[0]);}
diff --git a/src/ExprNode.cc b/src/ExprNode.cc
index de3cd447bf171b7905c38605eebc15fde93e2db7..a0d66216c05289278620aee6a2f1369b3304381a 100644
--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
@@ -782,7 +782,7 @@ VariableNode::prepareForDerivation()
       exit(EXIT_FAILURE);
     case SymbolType::excludedVariable:
       cerr << "VariableNode::prepareForDerivation: impossible case: "
-           << "You are trying to derive a variable that has been excluded via include_eqs/exclude_eqs: "
+           << "You are trying to derive a variable that has been excluded via model_remove/include_eqs/exclude_eqs: "
            << datatree.symbol_table.getName(symb_id) << endl;
       exit(EXIT_FAILURE);
     }
diff --git a/src/ModelTree.cc b/src/ModelTree.cc
index 3a7f98ba837a81d04a06745d5b60621be6d0b499..bbbde0c2573b289dde98422d19c63d954bb0c7c3 100644
--- a/src/ModelTree.cc
+++ b/src/ModelTree.cc
@@ -1545,87 +1545,6 @@ ModelTree::addEquation(expr_t eq, int lineno)
   equations_lineno.push_back(lineno);
 }
 
-vector<int>
-ModelTree::includeExcludeEquations(set<pair<string, string>> &eqs, bool exclude_eqs,
-                                   vector<BinaryOpNode *> &equations, vector<int> &equations_lineno,
-                                   EquationTags &equation_tags, bool static_equations) const
-{
-  vector<int> excluded_vars;
-  if (equations.empty())
-    return excluded_vars;
-
-  // Get equation numbers of tags
-  set<int> tag_eqns;
-  for (auto it = eqs.begin(); it != eqs.end();)
-    if (auto tmp = equation_tags.getEqnsByTag(it->first, it->second);
-        !tmp.empty())
-      {
-        tag_eqns.insert(tmp.begin(), tmp.end());
-        it = eqs.erase(it);
-      }
-    else
-      ++it;
-
-  set<int> eqns;
-  if (exclude_eqs)
-    eqns = tag_eqns;
-  else
-    for (size_t i = 0; i < equations.size(); i++)
-      if (tag_eqns.find(i) == tag_eqns.end())
-        eqns.insert(i);
-
-  // remove from equations, equations_lineno, equation_tags
-  vector<BinaryOpNode *> new_eqns;
-  vector<int> new_equations_lineno;
-  map<int, int> old_eqn_num_2_new;
-  for (size_t i = 0; i < equations.size(); i++)
-    if (eqns.find(i) != eqns.end())
-      {
-        if (auto tmp = equation_tags.getTagValueByEqnAndKey(i, "endogenous"); !tmp.empty())
-          excluded_vars.push_back(symbol_table.getID(tmp));
-        else
-          {
-            set<pair<int, int>> result;
-            equations[i]->arg1->collectDynamicVariables(SymbolType::endogenous, result);
-            if (result.size() == 1)
-              excluded_vars.push_back(result.begin()->first);
-            else
-              {
-                cerr << "ERROR: Equation " << i
-                     << " has been excluded but does not have a single variable on LHS or `endogenous` tag" << endl;
-                exit(EXIT_FAILURE);
-              }
-          }
-      }
-    else
-      {
-        new_eqns.emplace_back(equations[i]);
-        old_eqn_num_2_new[i] = new_eqns.size() - 1;
-        new_equations_lineno.emplace_back(equations_lineno[i]);
-      }
-  int n_excl = equations.size() - new_eqns.size();
-
-  equations = new_eqns;
-  equations_lineno = new_equations_lineno;
-
-  equation_tags.erase(eqns, old_eqn_num_2_new);
-
-  if (!static_equations)
-    for (size_t i = 0; i < excluded_vars.size(); i++)
-      for (size_t j = i+1; j < excluded_vars.size(); j++)
-        if (excluded_vars[i] == excluded_vars[j])
-          {
-            cerr << "Error: Variable " << symbol_table.getName(i) << " was excluded twice"
-                 << " via in/exclude_eqs option" << endl;
-            exit(EXIT_FAILURE);
-          }
-
-  cout << "Excluded " << n_excl << (static_equations ? " static " : " dynamic ")
-       << "equation" << (n_excl > 1 ? "s" : "") << " via in/exclude_eqs option" << endl;
-
-  return excluded_vars;
-}
-
 void
 ModelTree::findConstantEquationsWithoutMcpTag(map<VariableNode *, NumConstNode *> &subst_table) const
 {
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index 4596df6600cdd7e99d222a3c654b581b5353db95..373273e6b11708fc339c5934f5d4078e2b8658de 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -331,10 +331,6 @@ protected:
   void printBlockDecomposition() const;
   //! Determine for each block if it is linear or not
   void determineLinearBlocks();
-  //! Remove equations specified by exclude_eqs
-  vector<int> includeExcludeEquations(set<pair<string, string>> &eqs, bool exclude_eqs,
-                                      vector<BinaryOpNode *> &equations, vector<int> &equations_lineno,
-                                      EquationTags &equation_tags, bool static_equations) const;
 
   //! Return the type of equation belonging to the block
   EquationType
diff --git a/src/ParsingDriver.cc b/src/ParsingDriver.cc
index 23accb3389b0843027827aa4fb0ec67e41474558..a1363104175f4b9f8139a7489ee637e5923d3570 100644
--- a/src/ParsingDriver.cc
+++ b/src/ParsingDriver.cc
@@ -318,6 +318,8 @@ ParsingDriver::add_model_variable(const string &name)
   try
     {
       symb_id = mod_file->symbol_table.getID(name);
+      if (mod_file->symbol_table.getType(symb_id) == SymbolType::excludedVariable)
+        error("Variable '" + name + "' can no longer be used since it has been excluded by a previous 'model_remove' statement");
     }
   catch (SymbolTable::UnknownSymbolNameException &e)
     {
@@ -3529,3 +3531,16 @@ ParsingDriver::isSymbolIdentifier(const string &str)
       return false;
   return true;
 }
+
+void
+ParsingDriver::model_remove(const vector<pair<string, string>> &listed_eqs_by_tags)
+{
+  mod_file->dynamic_model.removeEquations(listed_eqs_by_tags, true, true);
+}
+
+void
+ParsingDriver::begin_model_replace(const vector<pair<string, string>> &listed_eqs_by_tags)
+{
+  mod_file->dynamic_model.removeEquations(listed_eqs_by_tags, true, false);
+  set_current_data_tree(&mod_file->dynamic_model);
+}
diff --git a/src/ParsingDriver.hh b/src/ParsingDriver.hh
index 35ba70ab0ae3ad9bb40786fa479050b5e9766587..61a55657d642c74a267cef1310d4d565e3805e7e 100644
--- a/src/ParsingDriver.hh
+++ b/src/ParsingDriver.hh
@@ -430,7 +430,7 @@ public:
   void add_epilogue_equal(const string &varname, expr_t expr);
   //! Begin a model block
   void begin_model();
-  //! End a model block, printing errors that were encountered in parsing
+  //! End a model or model_replace block, printing errors that were encountered in parsing
   void end_model();
   //! Writes a shocks statement
   void end_shocks(bool overwrite);
@@ -889,6 +889,10 @@ public:
   void begin_occbin_constraints();
   //! Add an occbin_constraints block
   void end_occbin_constraints(const vector<tuple<string, BinaryOpNode *, BinaryOpNode *, expr_t, expr_t>> &constraints);
+  // Process a model_remove statement
+  void model_remove(const vector<pair<string, string>> &listed_eqs_by_tags);
+  // Begin a model_replace statement
+  void begin_model_replace(const vector<pair<string, string>> &listed_eqs_by_tags);
   // Equivalent of MATLAB’s strsplit. Returns an empty vector given an empty string.
   static vector<string> strsplit(const string &str, char delim);
   // Returns true iff the string is a legal symbol identifier (see NAME token in lexer)