diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 54d823ade62e74df2245779e3984c633cb368d37..54a63e123a56e4e89ab2ec4abff5f44c3889eb6c 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -3353,10 +3353,14 @@ DynamicModel::fillVarModelTable() const
       for (const auto &eqtag : it.second)
         {
           set<pair<int, int>> lhs_set, lhs_tmp_set, rhs_set;
-          int eqn = equation_tags.getEqnByTag("name", eqtag);
-          if (eqn == -1)
+          int eqn;
+          try
+            {
+              eqn = equation_tags.getEqnByTag("name", eqtag);
+            }
+          catch (EquationTags::TagNotFoundException &e)
             {
-              cerr << "ERROR: equation tag '" << eqtag << "' not found" << endl;
+              cerr << "ERROR: no equation is named '" << eqtag << "'" << endl;
               exit(EXIT_FAILURE);
             }
 
@@ -3408,7 +3412,7 @@ DynamicModel::fillVarModelTable() const
   var_model_table.setLhsExprT(lhs_expr_tr);
 
   // Fill AR Matrix
-  var_model_table.setAR(fillAutoregressiveMatrix(true));
+  var_model_table.setAR(computeAutoregressiveMatrices(true));
 }
 
 void
@@ -3482,7 +3486,7 @@ DynamicModel::fillVarModelTableFromOrigModel() const
 }
 
 map<string, map<tuple<int, int, int>, expr_t>>
-DynamicModel::fillAutoregressiveMatrix(bool is_var) const
+DynamicModel::computeAutoregressiveMatrices(bool is_var) const
 {
   map<string, map<tuple<int, int, int>, expr_t>> ARr;
   auto eqnums = is_var ?
@@ -3516,10 +3520,14 @@ DynamicModel::fillTrendComponentModelTable() const
       vector<int> trend_eqnumber;
       for (const auto &eqtag : it.second)
         {
-          int eqn = equation_tags.getEqnByTag("name", eqtag);
-          if (eqn == -1)
+          int eqn;
+          try
             {
-              cerr << "ERROR: trend equation tag '" << eqtag << "' not found" << endl;
+              eqn = equation_tags.getEqnByTag("name", eqtag);
+            }
+          catch (EquationTags::TagNotFoundException &e)
+            {
+              cerr << "ERROR: no equation is named '" << eqtag << "'" << endl;
               exit(EXIT_FAILURE);
             }
           trend_eqnumber.push_back(eqn);
@@ -3536,10 +3544,14 @@ DynamicModel::fillTrendComponentModelTable() const
       for (const auto &eqtag : it.second)
         {
           set<pair<int, int>> lhs_set, lhs_tmp_set, rhs_set;
-          int eqn = equation_tags.getEqnByTag("name", eqtag);
-          if (eqn == -1)
+          int eqn;
+          try
+            {
+              eqn = equation_tags.getEqnByTag("name", eqtag);
+            }
+          catch (EquationTags::TagNotFoundException &e)
             {
-              cerr << "ERROR: equation tag '" << eqtag << "' not found" << endl;
+              cerr << "ERROR: no equation is named '" << eqtag << "'" << endl;
               exit(EXIT_FAILURE);
             }
 
@@ -3590,7 +3602,7 @@ DynamicModel::fillTrendComponentModelTable() const
 }
 
 pair<map<string, map<tuple<int, int, int>, expr_t>>, map<string, map<tuple<int, int, int>, expr_t>>>
-DynamicModel::fillErrorComponentMatrix(const ExprNode::subst_table_t &diff_subst_table) const
+DynamicModel::computeErrorComponentMatrices(const ExprNode::subst_table_t &diff_subst_table) const
 {
   map<string, map<tuple<int, int, int>, expr_t>> A0r, A0starr;
 
@@ -3694,11 +3706,11 @@ DynamicModel::fillTrendComponentModelTableFromOrigModel() const
 }
 
 void
-DynamicModel::fillTrendComponentmodelTableAREC(const ExprNode::subst_table_t &diff_subst_table) const
+DynamicModel::fillTrendComponentModelTableAREC(const ExprNode::subst_table_t &diff_subst_table) const
 {
-  auto ARr = fillAutoregressiveMatrix(false);
+  auto ARr = computeAutoregressiveMatrices(false);
   trend_component_model_table.setAR(ARr);
-  auto [A0r, A0starr] = fillErrorComponentMatrix(diff_subst_table);
+  auto [A0r, A0starr] = computeErrorComponentMatrices(diff_subst_table);
   trend_component_model_table.setA0(A0r, A0starr);
 }
 
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index 2031c596fd2d24bf648608b69f2c6e2a21d2baa0..ffe5e1a2556c5778bf2d8b33de4a827932309c03 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -271,6 +271,13 @@ private:
     str.erase(str.find_last_not_of("\t\n\v\f\r ") + 1);
   }
 
+  //! Compute autoregressive matrices of VAR or trend component models
+  map<string, map<tuple<int, int, int>, expr_t>> computeAutoregressiveMatrices(bool is_var) const;
+
+  //! Compute error component matrices of trend component_models
+  /*! Returns a pair (A0r, A0starr) */
+  pair<map<string, map<tuple<int, int, int>, expr_t>>, map<string, map<tuple<int, int, int>, expr_t>>> computeErrorComponentMatrices(const ExprNode::subst_table_t &diff_subst_table) const;
+
 public:
   DynamicModel(SymbolTable &symbol_table_arg,
                NumericalConstants &num_constants_arg,
@@ -345,20 +352,19 @@ public:
     return nonzero_hessian_eqs;
   }
 
-  //! Fill Autoregressive Matrix for var_model
-  map<string, map<tuple<int, int, int>, expr_t>> fillAutoregressiveMatrix(bool is_var) const;
-
-  //! Fill Error Component Matrix for trend_component_model
-  /*! Returns a pair (A0r, A0starr) */
-  pair<map<string, map<tuple<int, int, int>, expr_t>>, map<string, map<tuple<int, int, int>, expr_t>>> fillErrorComponentMatrix(const ExprNode::subst_table_t &diff_subst_table) const;
-
-  //! Fill the Trend Component Model Table
+  //! Fill the trend component model table with information available from the transformed model
   void fillTrendComponentModelTable() const;
+  //! Fill the trend component model table with information available from the original model
   void fillTrendComponentModelTableFromOrigModel() const;
-  void fillTrendComponentmodelTableAREC(const ExprNode::subst_table_t &diff_subst_table) const;
+  /* Fill the trend component model table with information about AR/EC
+     components, available from the transformed model. Needs to be called after
+     fillTrendComponentModelTableFromOrigModel() has been called on the
+     original model */
+  void fillTrendComponentModelTableAREC(const ExprNode::subst_table_t &diff_subst_table) const;
 
-  //! Fill the Var Model Table
+  //! Fill the VAR model table with information available from the transformed model
   void fillVarModelTable() const;
+  //! Fill the VAR model table with information available from the original model
   void fillVarModelTableFromOrigModel() const;
 
   //! Update the rhs references in the var model and trend component tables
diff --git a/src/EquationTags.cc b/src/EquationTags.cc
index a6693e4a0be63099ea5025405be71e3c84bba301..2d271375ddcfd7a273ea16a4e7a32a8f6318d6ab 100644
--- a/src/EquationTags.cc
+++ b/src/EquationTags.cc
@@ -47,7 +47,7 @@ EquationTags::getEqnByTag(const string &key, const string &value) const
   for (const auto & [eqn, tags] : eqn_tags)
     if (auto tmp = tags.find(key); tmp != tags.end() && tmp->second == value)
       return eqn;
-  return -1;
+  throw TagNotFoundException(key, value);
 }
 
 void
diff --git a/src/EquationTags.hh b/src/EquationTags.hh
index 57398bb5142c75b2ff132b2215831833d4299dbe..78d1db71c03c424f2b85e3803db4a9c00010e82b 100644
--- a/src/EquationTags.hh
+++ b/src/EquationTags.hh
@@ -31,6 +31,16 @@ class EquationTags
 private:
   map<int, map<string, string>> eqn_tags;
 public:
+  class TagNotFoundException
+  {
+  public:
+    const string key, value;
+    explicit TagNotFoundException(string key_arg, string value_arg)
+      : key{move(key_arg)}, value{move(value_arg)}
+    {
+    }
+  };
+
   // Add multiple equation tags for the given equation
   inline void add(int eqn, map<string, string> tags)
   {
@@ -89,7 +99,15 @@ public:
   //! Returns true if equation tag with key and value exists
   inline bool exists(const string &key, const string &value) const
   {
-    return getEqnByTag(key, value) >= 0;
+    try
+      {
+        getEqnByTag(key, value);
+      }
+    catch (TagNotFoundException &e)
+      {
+        return false;
+      }
+    return true;
   }
 
   inline bool exists(const int eqn) const
diff --git a/src/ModFile.cc b/src/ModFile.cc
index 00e0d2dfad7cc7d85fd2a9fbd9730f19061d40ba..78b6f88289d982f2209da2c83f6c8ac930976354 100644
--- a/src/ModFile.cc
+++ b/src/ModFile.cc
@@ -442,10 +442,10 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, bool
   // Create auxiliary variable and equations for Diff operators
   auto [diff_nodes, diff_subst_table] = dynamic_model.substituteDiff(pac_growth);
 
-  // Fill Trend Component Model Table
+  // Fill trend component and VAR model tables
   dynamic_model.fillTrendComponentModelTable();
   original_model.fillTrendComponentModelTableFromOrigModel();
-  dynamic_model.fillTrendComponentmodelTableAREC(diff_subst_table);
+  dynamic_model.fillTrendComponentModelTableAREC(diff_subst_table);
   dynamic_model.fillVarModelTable();
   original_model.fillVarModelTableFromOrigModel();
 
diff --git a/src/SubModel.cc b/src/SubModel.cc
index dfe0379b3d6a888003052b19416c63ef9cf29958..a41810a5351dab4df5968416ebbe31ab06aac805 100644
--- a/src/SubModel.cc
+++ b/src/SubModel.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2018-2019 Dynare Team
+ * Copyright © 2018-2021 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -118,13 +118,13 @@ TrendComponentModelTable::setA0(map<string, map<tuple<int, int, int>, expr_t>> A
   A0star = move(A0star_arg);
 }
 
-map<string, vector<string>>
+const map<string, vector<string>> &
 TrendComponentModelTable::getEqTags() const
 {
   return eqtags;
 }
 
-vector<string>
+const vector<string> &
 TrendComponentModelTable::getEqTags(const string &name_arg) const
 {
   checkModelName(name_arg);
@@ -142,80 +142,80 @@ TrendComponentModelTable::checkModelName(const string &name_arg) const
     }
 }
 
-vector<int>
+const vector<int> &
 TrendComponentModelTable::getNonTargetLhs(const string &name_arg) const
 {
   checkModelName(name_arg);
   return nontarget_lhs.find(name_arg)->second;
 }
 
-vector<int>
+const vector<int> &
 TrendComponentModelTable::getTargetLhs(const string &name_arg) const
 {
   checkModelName(name_arg);
   return target_lhs.find(name_arg)->second;
 }
 
-vector<int>
+const vector<int> &
 TrendComponentModelTable::getLhs(const string &name_arg) const
 {
   checkModelName(name_arg);
   return lhs.find(name_arg)->second;
 }
 
-vector<expr_t>
+const vector<expr_t> &
 TrendComponentModelTable::getLhsExprT(const string &name_arg) const
 {
   checkModelName(name_arg);
   return lhs_expr_t.find(name_arg)->second;
 }
 
-map<string, vector<string>>
+const map<string, vector<string>> &
 TrendComponentModelTable::getTargetEqTags() const
 {
   return target_eqtags;
 }
 
-map<string, vector<int>>
+const map<string, vector<int>> &
 TrendComponentModelTable::getEqNums() const
 {
   return eqnums;
 }
 
-map<string, vector<int>>
+const map<string, vector<int>> &
 TrendComponentModelTable::getTargetEqNums() const
 {
   return target_eqnums;
 }
 
-vector<int>
+const vector<int> &
 TrendComponentModelTable::getTargetEqNums(const string &name_arg) const
 {
   checkModelName(name_arg);
   return target_eqnums.find(name_arg)->second;
 }
 
-map<string, vector<int>>
+const map<string, vector<int>> &
 TrendComponentModelTable::getNonTargetEqNums() const
 {
   return nontarget_eqnums;
 }
 
-vector<int>
+const vector<int> &
 TrendComponentModelTable::getNonTargetEqNums(const string &name_arg) const
 {
   checkModelName(name_arg);
   return nontarget_eqnums.find(name_arg)->second;
 }
 
-vector<int>
+const vector<int> &
 TrendComponentModelTable::getEqNums(const string &name_arg) const
 {
   checkModelName(name_arg);
   return eqnums.find(name_arg)->second;
 }
 
-vector<int>
+const vector<int> &
 TrendComponentModelTable::getMaxLags(const string &name_arg) const
 {
   checkModelName(name_arg);
@@ -231,14 +231,14 @@ TrendComponentModelTable::getMaxLag(const string &name_arg) const
   return max_lag_int;
 }
 
-vector<bool>
+const vector<bool> &
 TrendComponentModelTable::getDiff(const string &name_arg) const
 {
   checkModelName(name_arg);
   return diff.find(name_arg)->second;
 }
 
-vector<int>
+const vector<int> &
 TrendComponentModelTable::getOrigDiffVar(const string &name_arg) const
 {
   checkModelName(name_arg);
@@ -306,12 +306,12 @@ TrendComponentModelTable::writeOutput(const string &basename, ostream &output) c
       for (const auto &it : rhs.at(name))
         {
           output << "M_.trend_component." << name << ".rhs.vars_at_eq{" << i << "}.var = [";
-          for (const auto &it1 : it)
-            output << symbol_table.getTypeSpecificID(it1.first) + 1 << " ";
+          for (auto [var, lag] : it)
+            output << symbol_table.getTypeSpecificID(var) + 1 << " ";
           output << "];" << endl
                  << "M_.trend_component." << name << ".rhs.vars_at_eq{" << i << "}.lag = [";
-          for (const auto &it1 : it)
-            output << it1.second << " ";
+          for (auto [var, lag] : it)
+            output << lag << " ";
           output << "];" << endl;
 
           i++;
@@ -350,42 +350,40 @@ TrendComponentModelTable::writeOutput(const string &basename, ostream &output) c
       ar_ec_output << "if strcmp(model_name, '" << name << "')" << endl
                    << "    % AR" << endl
                    << "    AR = zeros(" << nontarget_lhs_vec.size() << ", " << nontarget_lhs_vec.size() << ", " << getMaxLag(name) << ");" << endl;
-      for (const auto &it : AR.at(name))
+      for (const auto &[key, expr] : AR.at(name))
         {
-          auto [eqn, lag, lhs_symb_id] = it.first;
+          auto [eqn, lag, lhs_symb_id] = key;
           int colidx = static_cast<int>(distance(nontarget_lhs_vec.begin(), find(nontarget_lhs_vec.begin(), nontarget_lhs_vec.end(), lhs_symb_id)));
           ar_ec_output << "    AR(" << eqn + 1 << ", " << colidx + 1 << ", " << lag << ") = ";
-          it.second->writeOutput(ar_ec_output, ExprNodeOutputType::matlabDynamicModel);
+          expr->writeOutput(ar_ec_output, ExprNodeOutputType::matlabDynamicModel);
           ar_ec_output << ";" << endl;
         }
 
       int a0_lag = 0;
-      for (const auto &it : A0.at(name))
-        if (get<1>(it.first) > a0_lag)
-          a0_lag = get<1>(it.first);
+      for (const auto &[key, expr] : A0.at(name))
+        a0_lag = max(a0_lag, get<1>(key));
       ar_ec_output << endl
                    << "    % A0" << endl
                    << "    A0 = zeros(" << nontarget_lhs_vec.size() << ", " << nontarget_lhs_vec.size() << ", " << a0_lag << ");" << endl;
-      for (const auto &it : A0.at(name))
+      for (const auto &[key, expr] : A0.at(name))
         {
-          auto [eqn, lag, colidx] = it.first;
+          auto [eqn, lag, colidx] = key;
           ar_ec_output << "    A0(" << eqn + 1 << ", " << colidx + 1 << ", " << lag << ") = ";
-          it.second->writeOutput(ar_ec_output, ExprNodeOutputType::matlabDynamicModel);
+          expr->writeOutput(ar_ec_output, ExprNodeOutputType::matlabDynamicModel);
           ar_ec_output << ";" << endl;
         }
 
       int a0star_lag = 0;
-      for (const auto &it : A0star.at(name))
-        if (get<1>(it.first) > a0star_lag)
-          a0star_lag = get<1>(it.first);
+      for (const auto &[key, expr] : A0star.at(name))
+        a0star_lag = max(a0star_lag, get<1>(key));
       ar_ec_output << endl
                    << "    % A0star" << endl
                    << "    A0star = zeros(" << nontarget_lhs_vec.size() << ", " << target_lhs_vec.size() << ", " << a0star_lag << ");" << endl;
-      for (const auto &it : A0star.at(name))
+      for (const auto &[key, expr] : A0star.at(name))
         {
-          auto [eqn, lag, colidx] = it.first;
+          auto [eqn, lag, colidx] = key;
           ar_ec_output << "    A0star(" << eqn + 1 << ", " << colidx + 1 << ", " << lag << ") = ";
-          it.second->writeOutput(ar_ec_output, ExprNodeOutputType::matlabDynamicModel);
+          expr->writeOutput(ar_ec_output, ExprNodeOutputType::matlabDynamicModel);
           ar_ec_output << ";" << endl;
         }
 
@@ -444,7 +442,7 @@ VarModelTable::addVarModel(string name_arg, vector<string> eqtags_arg,
   names.insert(move(name_arg));
 }
 
-map<string, pair<SymbolList, int>>
+const map<string, pair<SymbolList, int>> &
 VarModelTable::getSymbolListAndOrder() const
 {
   return symbol_list_and_order;
@@ -471,11 +469,11 @@ VarModelTable::writeOutput(const string &basename, ostream &output) const
   for (const auto &name : names)
     {
       output << "M_.var." << name << ".model_name = '" << name << "';" << endl;
-      if (!symbol_list_and_order.at(name).first.empty())
+      if (auto &[symbol_list, order] = symbol_list_and_order.at(name);
+          !symbol_list.empty())
         {
-          symbol_list_and_order.at(name).first.writeOutput("M_.var." + name + ".var_list_", output);
-          output << "M_.var." << name << ".order = "
-                 << symbol_list_and_order.at(name).second << ";" << endl;
+          symbol_list.writeOutput("M_.var." + name + ".var_list_", output);
+          output << "M_.var." << name << ".order = " << order << ";" << endl;
         }
       output << "M_.var." << name << ".eqtags = {";
       for (const auto &it : eqtags.at(name))
@@ -506,12 +504,12 @@ VarModelTable::writeOutput(const string &basename, ostream &output) const
       for (const auto &it : rhs.at(name))
         {
           output << "M_.var." << name << ".rhs.vars_at_eq{" << i << "}.var = [";
-          for (const auto &it1 : it)
-            output << symbol_table.getTypeSpecificID(it1.first) + 1 << " ";
+          for (auto [var, lag] : it)
+            output << symbol_table.getTypeSpecificID(var) + 1 << " ";
           output << "];" << endl
                  << "M_.var." << name << ".rhs.vars_at_eq{" << i << "}.lag = [";
-          for (const auto &it1 : it)
-            output << it1.second << " ";
+          for (auto [var, lag] : it)
+            output << lag << " ";
           output << "];" << endl;
 
           i++;
@@ -520,12 +518,12 @@ VarModelTable::writeOutput(const string &basename, ostream &output) const
       vector<int> lhs = getLhsOrigIds(name);
       ar_output << "if strcmp(model_name, '" << name << "')" << endl
                 << "    ar = zeros(" << lhs.size() << ", " << lhs.size() << ", " << getMaxLag(name) << ");" << endl;
-      for (const auto &it : AR.at(name))
+      for (const auto &[key, expr] : AR.at(name))
         {
-          auto [eqn, lag, lhs_symb_id] = it.first;
+          auto [eqn, lag, lhs_symb_id] = key;
           int colidx = static_cast<int>(distance(lhs.begin(), find(lhs.begin(), lhs.end(), lhs_symb_id)));
           ar_output << "    ar(" << eqn + 1 << ", " << colidx + 1 << ", " << lag << ") = ";
-          it.second->writeOutput(ar_output, ExprNodeOutputType::matlabDynamicModel);
+          expr->writeOutput(ar_output, ExprNodeOutputType::matlabDynamicModel);
           ar_output << ";" << endl;
         }
       ar_output << "    return" << endl
@@ -562,13 +560,13 @@ VarModelTable::writeJsonOutput(ostream &output) const
     }
 }
 
-map<string, vector<string>>
+const map<string, vector<string>> &
 VarModelTable::getEqTags() const
 {
   return eqtags;
 }
 
-vector<string>
+const vector<string> &
 VarModelTable::getEqTags(const string &name_arg) const
 {
   checkModelName(name_arg);
@@ -630,20 +628,20 @@ VarModelTable::setLhsExprT(map<string, vector<expr_t>> lhs_expr_t_arg)
   lhs_expr_t = move(lhs_expr_t_arg);
 }
 
-map<string, vector<int>>
+const map<string, vector<int>> &
 VarModelTable::getEqNums() const
 {
   return eqnums;
 }
 
-vector<bool>
+const vector<bool> &
 VarModelTable::getDiff(const string &name_arg) const
 {
   checkModelName(name_arg);
   return diff.find(name_arg)->second;
 }
 
-vector<int>
+const vector<int> &
 VarModelTable::getEqNums(const string &name_arg) const
 {
   checkModelName(name_arg);
@@ -674,7 +672,7 @@ VarModelTable::setAR(map<string, map<tuple<int, int, int>, expr_t>> AR_arg)
   AR = move(AR_arg);
 }
 
-vector<int>
+const vector<int> &
 VarModelTable::getMaxLags(const string &name_arg) const
 {
   checkModelName(name_arg);
@@ -690,28 +688,28 @@ VarModelTable::getMaxLag(const string &name_arg) const
   return max_lag_int;
 }
 
-vector<int>
+const vector<int> &
 VarModelTable::getLhs(const string &name_arg) const
 {
   checkModelName(name_arg);
   return lhs.find(name_arg)->second;
 }
 
-vector<int>
+const vector<int> &
 VarModelTable::getLhsOrigIds(const string &name_arg) const
 {
   checkModelName(name_arg);
   return lhs_orig_symb_ids.find(name_arg)->second;
 }
 
-vector<set<pair<int, int>>>
+const vector<set<pair<int, int>>> &
 VarModelTable::getRhs(const string &name_arg) const
 {
   checkModelName(name_arg);
   return rhs.find(name_arg)->second;
 }
 
-vector<expr_t>
+const vector<expr_t> &
 VarModelTable::getLhsExprT(const string &name_arg) const
 {
   checkModelName(name_arg);
diff --git a/src/SubModel.hh b/src/SubModel.hh
index 415ae14ea579657cf5c881b5c749f6226566cdd8..b30a0c2a00aa66f7d4af038d939ba00831e50b83 100644
--- a/src/SubModel.hh
+++ b/src/SubModel.hh
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2018-2019 Dynare Team
+ * Copyright © 2018-2021 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -46,8 +46,8 @@ private:
   map<string, vector<bool>> diff;
   map<string, vector<expr_t>> lhs_expr_t;
   map<string, vector<int>> target_vars;
-  map<string, map<tuple<int, int, int>, expr_t>> AR; // AR: name -> (eqn, lag, lhs_symb_id) -> expr_t
-  map<string, map<tuple<int, int, int>, expr_t>> A0, A0star; // EC: name -> (eqn, lag, col) -> expr_t
+  map<string, map<tuple<int, int, int>, expr_t>> AR; // name -> (eqn, lag, lhs_symb_id) -> expr_t
+  map<string, map<tuple<int, int, int>, expr_t>> A0, A0star; // name -> (eqn, lag, col) -> expr_t
 public:
   explicit TrendComponentModelTable(SymbolTable &symbol_table_arg);
 
@@ -58,23 +58,23 @@ public:
   inline bool isExistingTrendComponentModelName(const string &name_arg) const;
   inline bool empty() const;
 
-  map<string, vector<string>> getEqTags() const;
-  vector<string> getEqTags(const string &name_arg) const;
-  map<string, vector<string>> getTargetEqTags() const;
-  map<string, vector<int>> getEqNums() const;
-  map<string, vector<int>> getTargetEqNums() const;
-  vector<int> getTargetEqNums(const string &name_arg) const;
-  vector<int> getEqNums(const string &name_arg) const;
-  vector<int> getMaxLags(const string &name_arg) const;
+  const map<string, vector<string>> &getEqTags() const;
+  const vector<string> &getEqTags(const string &name_arg) const;
+  const map<string, vector<string>> &getTargetEqTags() const;
+  const map<string, vector<int>> &getEqNums() const;
+  const map<string, vector<int>> &getTargetEqNums() const;
+  const vector<int> &getTargetEqNums(const string &name_arg) const;
+  const vector<int> &getEqNums(const string &name_arg) const;
+  const vector<int> &getMaxLags(const string &name_arg) const;
   int getMaxLag(const string &name_arg) const;
-  vector<int> getLhs(const string &name_arg) const;
-  vector<expr_t> getLhsExprT(const string &name_arg) const;
-  vector<bool> getDiff(const string &name_arg) const;
-  vector<int> getOrigDiffVar(const string &name_arg) const;
-  map<string, vector<int>> getNonTargetEqNums() const;
-  vector<int> getNonTargetEqNums(const string &name_arg) const;
-  vector<int> getNonTargetLhs(const string &name_arg) const;
-  vector<int> getTargetLhs(const string &name_arg) const;
+  const vector<int> &getLhs(const string &name_arg) const;
+  const vector<expr_t> &getLhsExprT(const string &name_arg) const;
+  const vector<bool> &getDiff(const string &name_arg) const;
+  const vector<int> &getOrigDiffVar(const string &name_arg) const;
+  const map<string, vector<int>> &getNonTargetEqNums() const;
+  const vector<int> &getNonTargetEqNums(const string &name_arg) const;
+  const vector<int> &getNonTargetLhs(const string &name_arg) const;
+  const vector<int> &getTargetLhs(const string &name_arg) const;
 
   void setVals(map<string, vector<int>> eqnums_arg, map<string, vector<int>> target_eqnums_arg,
                map<string, vector<int>> lhs_arg,
@@ -119,10 +119,10 @@ private:
   map<string, pair<SymbolList, int>> symbol_list_and_order;
   map<string, vector<string>> eqtags;
   map<string, vector<int>> eqnums, max_lags, lhs, lhs_orig_symb_ids, orig_diff_var;
-  map<string, vector<set<pair<int, int>>>> rhs;
+  map<string, vector<set<pair<int, int>>>> rhs; // name -> for each equation: set of pairs (var, lag)
   map<string, vector<bool>> diff;
   map<string, vector<expr_t>> lhs_expr_t;
-  map<string, map<tuple<int, int, int>, expr_t>> AR; // AR: name -> (eqn, lag, lhs_symb_id) -> param_expr_t
+  map<string, map<tuple<int, int, int>, expr_t>> AR; // name -> (eqn, lag, lhs_symb_id) -> param_expr_t
 public:
   explicit VarModelTable(SymbolTable &symbol_table_arg);
 
@@ -133,18 +133,18 @@ public:
   inline bool isExistingVarModelName(const string &name_arg) const;
   inline bool empty() const;
 
-  map<string, vector<string>> getEqTags() const;
-  vector<string> getEqTags(const string &name_arg) const;
-  map<string, vector<int>> getEqNums() const;
-  vector<bool> getDiff(const string &name_arg) const;
-  vector<int> getEqNums(const string &name_arg) const;
-  vector<int> getMaxLags(const string &name_arg) const;
+  const map<string, vector<string>> &getEqTags() const;
+  const vector<string> &getEqTags(const string &name_arg) const;
+  const map<string, vector<int>> &getEqNums() const;
+  const vector<bool> &getDiff(const string &name_arg) const;
+  const vector<int> &getEqNums(const string &name_arg) const;
+  const vector<int> &getMaxLags(const string &name_arg) const;
   int getMaxLag(const string &name_arg) const;
-  vector<int> getLhs(const string &name_arg) const;
-  vector<int> getLhsOrigIds(const string &name_arg) const;
-  map<string, pair<SymbolList, int>> getSymbolListAndOrder() const;
-  vector<set<pair<int, int>>> getRhs(const string &name_arg) const;
-  vector<expr_t> getLhsExprT(const string &name_arg) const;
+  const vector<int> &getLhs(const string &name_arg) const;
+  const vector<int> &getLhsOrigIds(const string &name_arg) const;
+  const map<string, pair<SymbolList, int>> &getSymbolListAndOrder() const;
+  const vector<set<pair<int, int>>> &getRhs(const string &name_arg) const;
+  const vector<expr_t> &getLhsExprT(const string &name_arg) const;
 
   void setEqNums(map<string, vector<int>> eqnums_arg);
   void setLhs(map<string, vector<int>> lhs_arg);