diff --git a/src/ComputingTasks.cc b/src/ComputingTasks.cc
index f3eade18..9747055f 100644
--- a/src/ComputingTasks.cc
+++ b/src/ComputingTasks.cc
@@ -1047,17 +1047,16 @@ RamseyPolicyStatement::checkPass(ModFileStructure &mod_file_struct, WarningConso
 void
 RamseyPolicyStatement::checkRamseyPolicyList()
 {
-  for (vector<string>::const_iterator it = ramsey_policy_list.begin();
-       it != ramsey_policy_list.end(); it++)
+  for (const auto & it : ramsey_policy_list)
     {
-      if (!symbol_table.exists(*it))
+      if (!symbol_table.exists(it))
         {
-          cerr << "ERROR: ramsey_policy: " << *it << " was not declared." << endl;
+          cerr << "ERROR: ramsey_policy: " << it << " was not declared." << endl;
           exit(EXIT_FAILURE);
         }
-      if (symbol_table.getType(*it) != eEndogenous)
+      if (symbol_table.getType(it) != eEndogenous)
         {
-          cerr << "ERROR: ramsey_policy: " << *it << " is not endogenous." << endl;
+          cerr << "ERROR: ramsey_policy: " << it << " is not endogenous." << endl;
           exit(EXIT_FAILURE);
         }
     }
@@ -1449,18 +1448,17 @@ EstimatedParamsStatement::EstimatedParamsStatement(const vector<EstimationParams
 void
 EstimatedParamsStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings)
 {
-  for (vector<EstimationParams>::const_iterator it = estim_params_list.begin();
-       it != estim_params_list.end(); it++)
+  for (const auto & it : estim_params_list)
     {
-      if (it->name == "dsge_prior_weight")
+      if (it.name == "dsge_prior_weight")
         mod_file_struct.dsge_prior_weight_in_estimated_params = true;
 
       // Handle case of degenerate beta prior
-      if (it->prior == eBeta)
+      if (it.prior == eBeta)
         try
           {
-            if (it->mean->eval(eval_context_t()) == 0.5
-                && it->std->eval(eval_context_t()) == 0.5)
+            if (it.mean->eval(eval_context_t()) == 0.5
+                && it.std->eval(eval_context_t()) == 0.5)
               {
                 cerr << "ERROR: The prior density is not defined for the beta distribution when the mean = standard deviation = 0.5." << endl;
                 exit(EXIT_FAILURE);
@@ -1475,39 +1473,37 @@ EstimatedParamsStatement::checkPass(ModFileStructure &mod_file_struct, WarningCo
   // Check that no parameter/endogenous is declared twice in the block
   set<string> already_declared;
   set<pair<string, string> > already_declared_corr;
-  for (vector<EstimationParams>::const_iterator it = estim_params_list.begin();
-       it != estim_params_list.end(); it++)
+  for (const auto & it : estim_params_list)
     {
-      if (it->type == 3) // Correlation
+      if (it.type == 3) // Correlation
         {
           // Use lexical ordering for the pair of symbols
-          pair<string, string> x = it->name < it->name2 ? make_pair(it->name, it->name2) : make_pair(it->name2, it->name);
+          pair<string, string> x = it.name < it.name2 ? make_pair(it.name, it.name2) : make_pair(it.name2, it.name);
 
           if (already_declared_corr.find(x) == already_declared_corr.end())
             already_declared_corr.insert(x);
           else
             {
-              cerr << "ERROR: in `estimated_params' block, the correlation between " << it->name << " and " << it->name2 << " is declared twice." << endl;
+              cerr << "ERROR: in `estimated_params' block, the correlation between " << it.name << " and " << it.name2 << " is declared twice." << endl;
               exit(EXIT_FAILURE);
             }
         }
       else
         {
-          if (already_declared.find(it->name) == already_declared.end())
-            already_declared.insert(it->name);
+          if (already_declared.find(it.name) == already_declared.end())
+            already_declared.insert(it.name);
           else
             {
-              cerr << "ERROR: in `estimated_params' block, the symbol " << it->name << " is declared twice." << endl;
+              cerr << "ERROR: in `estimated_params' block, the symbol " << it.name << " is declared twice." << endl;
               exit(EXIT_FAILURE);
             }
         }
     }
 
   // Fill in mod_file_struct.estimated_parameters (related to #469)
-  for (vector<EstimationParams>::const_iterator it = estim_params_list.begin();
-       it != estim_params_list.end(); it++)
-    if (it->type == 2 && it->name != "dsge_prior_weight")
-      mod_file_struct.estimated_parameters.insert(symbol_table.getID(it->name));
+  for (const auto & it : estim_params_list)
+    if (it.type == 2 && it.name != "dsge_prior_weight")
+      mod_file_struct.estimated_parameters.insert(symbol_table.getID(it.name));
 }
 
 void
@@ -1519,12 +1515,12 @@ EstimatedParamsStatement::writeOutput(ostream &output, const string &basename, b
          << "estim_params_.corrn = [];" << endl
          << "estim_params_.param_vals = [];" << endl;
 
-  for (vector<EstimationParams>::const_iterator it = estim_params_list.begin(); it != estim_params_list.end(); it++)
+  for (const auto & it : estim_params_list)
     {
-      int symb_id = symbol_table.getTypeSpecificID(it->name) + 1;
-      SymbolType symb_type = symbol_table.getType(it->name);
+      int symb_id = symbol_table.getTypeSpecificID(it.name) + 1;
+      SymbolType symb_type = symbol_table.getType(it.name);
 
-      switch (it->type)
+      switch (it.type)
         {
         case 1:
           if (symb_type == eExogenous)
@@ -1542,26 +1538,26 @@ EstimatedParamsStatement::writeOutput(ostream &output, const string &basename, b
             output << "estim_params_.corrx = [estim_params_.corrx; ";
           else if (symb_type == eEndogenous)
             output << "estim_params_.corrn = [estim_params_.corrn; ";
-          output << symb_id << " " << symbol_table.getTypeSpecificID(it->name2)+1;
+          output << symb_id << " " << symbol_table.getTypeSpecificID(it.name2)+1;
           break;
         }
       output << ", ";
-      it->init_val->writeOutput(output);
+      it.init_val->writeOutput(output);
       output << ", ";
-      it->low_bound->writeOutput(output);
+      it.low_bound->writeOutput(output);
       output << ", ";
-      it->up_bound->writeOutput(output);
+      it.up_bound->writeOutput(output);
       output << ", "
-             << it->prior << ", ";
-      it->mean->writeOutput(output);
+             << it.prior << ", ";
+      it.mean->writeOutput(output);
       output << ", ";
-      it->std->writeOutput(output);
+      it.std->writeOutput(output);
       output << ", ";
-      it->p3->writeOutput(output);
+      it.p3->writeOutput(output);
       output << ", ";
-      it->p4->writeOutput(output);
+      it.p4->writeOutput(output);
       output << ", ";
-      it->jscale->writeOutput(output);
+      it.jscale->writeOutput(output);
       output << " ];" << endl;
     }
 }
@@ -1636,32 +1632,32 @@ EstimatedParamsInitStatement::writeOutput(ostream &output, const string &basenam
   if (use_calibration)
     output << "options_.use_calibration_initialization = 1;" << endl;
 
-  for (vector<EstimationParams>::const_iterator it = estim_params_list.begin(); it != estim_params_list.end(); it++)
+  for (const auto & it : estim_params_list)
     {
-      int symb_id = symbol_table.getTypeSpecificID(it->name) + 1;
-      SymbolType symb_type = symbol_table.getType(it->name);
+      int symb_id = symbol_table.getTypeSpecificID(it.name) + 1;
+      SymbolType symb_type = symbol_table.getType(it.name);
 
-      if (it->type < 3)
+      if (it.type < 3)
         {
           if (symb_type == eExogenous)
             {
               output << "tmp1 = find(estim_params_.var_exo(:,1)==" << symb_id << ");" << endl;
               output << "estim_params_.var_exo(tmp1,2) = ";
-              it->init_val->writeOutput(output);
+              it.init_val->writeOutput(output);
               output << ";" << endl;
             }
           else if (symb_type == eEndogenous)
             {
               output << "tmp1 = find(estim_params_.var_endo(:,1)==" << symb_id << ");" << endl;
               output << "estim_params_.var_endo(tmp1,2) = ";
-              it->init_val->writeOutput(output);
+              it.init_val->writeOutput(output);
               output << ";" << endl;
             }
           else if (symb_type == eParameter)
             {
               output << "tmp1 = find(estim_params_.param_vals(:,1)==" << symb_id << ");" << endl;
               output << "estim_params_.param_vals(tmp1,2) = ";
-              it->init_val->writeOutput(output);
+              it.init_val->writeOutput(output);
               output << ";" << endl;
             }
         }
@@ -1669,18 +1665,18 @@ EstimatedParamsInitStatement::writeOutput(ostream &output, const string &basenam
         {
           if (symb_type == eExogenous)
             {
-              output << "tmp1 = find((estim_params_.corrx(:,1)==" << symb_id << " & estim_params_.corrx(:,2)==" << symbol_table.getTypeSpecificID(it->name2)+1 << ") | "
-                     <<             "(estim_params_.corrx(:,2)==" << symb_id << " & estim_params_.corrx(:,1)==" << symbol_table.getTypeSpecificID(it->name2)+1 << "));" << endl;
+              output << "tmp1 = find((estim_params_.corrx(:,1)==" << symb_id << " & estim_params_.corrx(:,2)==" << symbol_table.getTypeSpecificID(it.name2)+1 << ") | "
+                     <<             "(estim_params_.corrx(:,2)==" << symb_id << " & estim_params_.corrx(:,1)==" << symbol_table.getTypeSpecificID(it.name2)+1 << "));" << endl;
               output << "estim_params_.corrx(tmp1,3) = ";
-              it->init_val->writeOutput(output);
+              it.init_val->writeOutput(output);
               output << ";" << endl;
             }
           else if (symb_type == eEndogenous)
             {
-              output << "tmp1 = find((estim_params_.corrn(:,1)==" << symb_id << " & estim_params_.corrn(:,2)==" << symbol_table.getTypeSpecificID(it->name2)+1 << ") | "
-                     <<             "(estim_params_.corrn(:,2)==" << symb_id << " & estim_params_.corrn(:,1)==" << symbol_table.getTypeSpecificID(it->name2)+1 << "));" << endl;
+              output << "tmp1 = find((estim_params_.corrn(:,1)==" << symb_id << " & estim_params_.corrn(:,2)==" << symbol_table.getTypeSpecificID(it.name2)+1 << ") | "
+                     <<             "(estim_params_.corrn(:,2)==" << symb_id << " & estim_params_.corrn(:,1)==" << symbol_table.getTypeSpecificID(it.name2)+1 << "));" << endl;
               output << "estim_params_.corrn(tmp1,3) = ";
-              it->init_val->writeOutput(output);
+              it.init_val->writeOutput(output);
               output << ";" << endl;
             }
         }
@@ -1732,23 +1728,23 @@ EstimatedParamsBoundsStatement::EstimatedParamsBoundsStatement(const vector<Esti
 void
 EstimatedParamsBoundsStatement::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const
 {
-  for (vector<EstimationParams>::const_iterator it = estim_params_list.begin(); it != estim_params_list.end(); it++)
+  for (const auto & it : estim_params_list)
     {
-      int symb_id = symbol_table.getTypeSpecificID(it->name) + 1;
-      SymbolType symb_type = symbol_table.getType(it->name);
+      int symb_id = symbol_table.getTypeSpecificID(it.name) + 1;
+      SymbolType symb_type = symbol_table.getType(it.name);
 
-      if (it->type < 3)
+      if (it.type < 3)
         {
           if (symb_type == eExogenous)
             {
               output << "tmp1 = find(estim_params_.var_exo(:,1)==" << symb_id << ");" << endl;
 
               output << "estim_params_.var_exo(tmp1,3) = ";
-              it->low_bound->writeOutput(output);
+              it.low_bound->writeOutput(output);
               output << ";" << endl;
 
               output << "estim_params_.var_exo(tmp1,4) = ";
-              it->up_bound->writeOutput(output);
+              it.up_bound->writeOutput(output);
               output << ";" << endl;
             }
           else if (symb_type == eEndogenous)
@@ -1756,11 +1752,11 @@ EstimatedParamsBoundsStatement::writeOutput(ostream &output, const string &basen
               output << "tmp1 = find(estim_params_.var_endo(:,1)==" << symb_id << ");" << endl;
 
               output << "estim_params_.var_endo(tmp1,3) = ";
-              it->low_bound->writeOutput(output);
+              it.low_bound->writeOutput(output);
               output << ";" << endl;
 
               output << "estim_params_.var_endo(tmp1,4) = ";
-              it->up_bound->writeOutput(output);
+              it.up_bound->writeOutput(output);
               output << ";" << endl;
             }
           else if (symb_type == eParameter)
@@ -1768,11 +1764,11 @@ EstimatedParamsBoundsStatement::writeOutput(ostream &output, const string &basen
               output << "tmp1 = find(estim_params_.param_vals(:,1)==" << symb_id << ");" << endl;
 
               output << "estim_params_.param_vals(tmp1,3) = ";
-              it->low_bound->writeOutput(output);
+              it.low_bound->writeOutput(output);
               output << ";" << endl;
 
               output << "estim_params_.param_vals(tmp1,4) = ";
-              it->up_bound->writeOutput(output);
+              it.up_bound->writeOutput(output);
               output << ";" << endl;
             }
         }
@@ -1780,28 +1776,28 @@ EstimatedParamsBoundsStatement::writeOutput(ostream &output, const string &basen
         {
           if (symb_type == eExogenous)
             {
-              output << "tmp1 = find((estim_params_.corrx(:,1)==" << symb_id << " & estim_params_.corrx(:,2)==" << symbol_table.getTypeSpecificID(it->name2)+1 << ") | "
-                     <<             "(estim_params_.corrx(:,2)==" << symb_id << " & estim_params_.corrx(:,1)==" << symbol_table.getTypeSpecificID(it->name2)+1 << "));" << endl;
+              output << "tmp1 = find((estim_params_.corrx(:,1)==" << symb_id << " & estim_params_.corrx(:,2)==" << symbol_table.getTypeSpecificID(it.name2)+1 << ") | "
+                     <<             "(estim_params_.corrx(:,2)==" << symb_id << " & estim_params_.corrx(:,1)==" << symbol_table.getTypeSpecificID(it.name2)+1 << "));" << endl;
 
               output << "estim_params_.corrx(tmp1,4) = ";
-              it->low_bound->writeOutput(output);
+              it.low_bound->writeOutput(output);
               output << ";" << endl;
 
               output << "estim_params_.corrx(tmp1,5) = ";
-              it->up_bound->writeOutput(output);
+              it.up_bound->writeOutput(output);
               output << ";" << endl;
             }
           else if (symb_type == eEndogenous)
             {
-              output << "tmp1 = find((estim_params_.corrn(:,1)==" << symb_id << " & estim_params_.corrn(:,2)==" << symbol_table.getTypeSpecificID(it->name2)+1 << ") | "
-                     <<             "(estim_params_.corrn(:,2)==" << symb_id << " & estim_params_.corrn(:,1)==" << symbol_table.getTypeSpecificID(it->name2)+1 << "));" << endl;
+              output << "tmp1 = find((estim_params_.corrn(:,1)==" << symb_id << " & estim_params_.corrn(:,2)==" << symbol_table.getTypeSpecificID(it.name2)+1 << ") | "
+                     <<             "(estim_params_.corrn(:,2)==" << symb_id << " & estim_params_.corrn(:,1)==" << symbol_table.getTypeSpecificID(it.name2)+1 << "));" << endl;
 
               output << "estim_params_.corrn(tmp1,4) = ";
-              it->low_bound->writeOutput(output);
+              it.low_bound->writeOutput(output);
               output << ";" << endl;
 
               output << "estim_params_.corrn(tmp1,5) = ";
-              it->up_bound->writeOutput(output);
+              it.up_bound->writeOutput(output);
               output << ";" << endl;
             }
         }
@@ -1852,18 +1848,18 @@ void
 ObservationTrendsStatement::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const
 {
   output << "options_.trend_coeff = {};" << endl;
-  for (trend_elements_t::const_iterator it = trend_elements.begin(); it != trend_elements.end(); it++)
+  for (const auto & trend_element : trend_elements)
     {
-      SymbolType type = symbol_table.getType(it->first);
+      SymbolType type = symbol_table.getType(trend_element.first);
       if (type == eEndogenous)
         {
-          output << "tmp1 = strmatch('" << it->first << "',options_.varobs,'exact');" << endl;
+          output << "tmp1 = strmatch('" << trend_element.first << "',options_.varobs,'exact');" << endl;
           output << "options_.trend_coeffs{tmp1} = '";
-          it->second->writeOutput(output);
+          trend_element.second->writeOutput(output);
           output << "';" << endl;
         }
       else
-        cerr << "Warning : Non-variable symbol used in observation_trends: " << it->first << endl;
+        cerr << "Warning : Non-variable symbol used in observation_trends: " << trend_element.first << endl;
     }
 }
 
@@ -1873,20 +1869,19 @@ ObservationTrendsStatement::writeJsonOutput(ostream &output) const
   output << "{\"statementName\": \"observation_trends\", "
          << "\"trends\" : {";
   bool printed = false;
-  for (trend_elements_t::const_iterator it = trend_elements.begin();
-       it != trend_elements.end(); it++)
+  for (const auto & trend_element : trend_elements)
     {
-      if (symbol_table.getType(it->first) == eEndogenous)
+      if (symbol_table.getType(trend_element.first) == eEndogenous)
         {
           if (printed)
             output << ", ";
-          output << "\"" << it->first << "\": \"";
-          it->second->writeJsonOutput(output, {}, {});
+          output << "\"" << trend_element.first << "\": \"";
+          trend_element.second->writeJsonOutput(output, {}, {});
           output << "\"" << endl;
           printed = true;
         }
       else
-        cerr << "Warning : Non-variable symbol used in observation_trends: " << it->first << endl;
+        cerr << "Warning : Non-variable symbol used in observation_trends: " << trend_element.first << endl;
     }
   output << "}"
          << "}";
@@ -1958,13 +1953,12 @@ OsrParamsBoundsStatement::writeOutput(ostream &output, const string &basename, b
 
   output << "M_.osr.param_bounds = [-inf(length(M_.osr.param_names), 1), inf(length(M_.osr.param_names), 1)];" << endl;
 
-  for (vector<OsrParams>::const_iterator it = osr_params_list.begin();
-       it != osr_params_list.end(); it++)
+  for (const auto & it : osr_params_list)
     {
-      output << "M_.osr.param_bounds(strcmp(M_.osr.param_names, '" << it->name << "'), :) = [";
-      it->low_bound->writeOutput(output);
+      output << "M_.osr.param_bounds(strcmp(M_.osr.param_names, '" << it.name << "'), :) = [";
+      it.low_bound->writeOutput(output);
       output << ", ";
-      it->up_bound->writeOutput(output);
+      it.up_bound->writeOutput(output);
       output << "];" << endl;
     }
 }
@@ -2069,11 +2063,10 @@ OptimWeightsStatement::writeOutput(ostream &output, const string &basename, bool
          << "M_.osr.variable_weights = sparse(M_.endo_nbr,M_.endo_nbr);" << endl
          << "M_.osr.variable_indices = [];" << endl << endl;
 
-  for (var_weights_t::const_iterator it = var_weights.begin();
-       it != var_weights.end(); it++)
+  for (const auto & var_weight : var_weights)
     {
-      const string &name = it->first;
-      const expr_t value = it->second;
+      const string &name = var_weight.first;
+      const expr_t value = var_weight.second;
       int id = symbol_table.getTypeSpecificID(name) + 1;
       output << "M_.osr.variable_weights(" << id << "," << id << ") = ";
       value->writeOutput(output);
@@ -2081,12 +2074,11 @@ OptimWeightsStatement::writeOutput(ostream &output, const string &basename, bool
       output << "M_.osr.variable_indices = [M_.osr.variable_indices; " << id << "];" << endl;
     }
 
-  for (covar_weights_t::const_iterator it = covar_weights.begin();
-       it != covar_weights.end(); it++)
+  for (const auto & covar_weight : covar_weights)
     {
-      const string &name1 = it->first.first;
-      const string &name2 = it->first.second;
-      const expr_t value = it->second;
+      const string &name1 = covar_weight.first.first;
+      const string &name2 = covar_weight.first.second;
+      const expr_t value = covar_weight.second;
       int id1 = symbol_table.getTypeSpecificID(name1) + 1;
       int id2 = symbol_table.getTypeSpecificID(name2) + 1;
       output << "M_.osr.variable_weights(" << id1 << "," << id2 << ") = ";
@@ -2198,11 +2190,10 @@ ModelComparisonStatement::writeOutput(ostream &output, const string &basename, b
   output << "ModelNames_ = {};" << endl;
   output << "ModelPriors_ = [];" << endl;
 
-  for (filename_list_t::const_iterator it = filename_list.begin();
-       it != filename_list.end(); it++)
+  for (const auto & it : filename_list)
     {
-      output << "ModelNames_ = { ModelNames_{:} '" << (*it).first << "'};" << endl;
-      output << "ModelPriors_ = [ ModelPriors_ ; " << (*it).second << "];" << endl;
+      output << "ModelNames_ = { ModelNames_{:} '" << it.first << "'};" << endl;
+      output << "ModelPriors_ = [ ModelPriors_ ; " << it.second << "];" << endl;
     }
   output << "oo_ = model_comparison(ModelNames_,ModelPriors_,oo_,options_,M_.fname);" << endl;
 }
@@ -2970,9 +2961,9 @@ int
 SvarIdentificationStatement::getMaxLag() const
 {
   int max_lag = 0;
-  for (svar_identification_restrictions_t::const_iterator it = restrictions.begin(); it != restrictions.end(); it++)
-    if (it->lag > max_lag)
-      max_lag = it->lag;
+  for (const auto & restriction : restrictions)
+    if (restriction.lag > max_lag)
+      max_lag = restriction.lag;
 
   return max_lag;
 }
@@ -3113,12 +3104,11 @@ MarkovSwitchingStatement::MarkovSwitchingStatement(const OptionsList &options_li
 
       vector<string> tokenizedRestrictions;
       split(tokenizedRestrictions, it_num->second, is_any_of("["), token_compress_on);
-      for (vector<string>::iterator it = tokenizedRestrictions.begin();
-           it != tokenizedRestrictions.end(); it++)
-        if (it->size() > 0)
+      for (auto & tokenizedRestriction : tokenizedRestrictions)
+        if (tokenizedRestriction.size() > 0)
           {
             vector<string> restriction;
-            split(restriction, *it, is_any_of("], "));
+            split(restriction, tokenizedRestriction, is_any_of("], "));
             for (vector<string>::iterator it1 = restriction.begin();
                  it1 != restriction.end();)
               if (it1->empty())
@@ -3165,7 +3155,7 @@ MarkovSwitchingStatement::MarkovSwitchingStatement(const OptionsList &options_li
               {
                 cerr << "ERROR: The first two arguments for a restriction must be integers "
                      << "specifying the regime and the last must be a double specifying the "
-                     << "transition probability. You wrote [" << *it << endl;
+                     << "transition probability. You wrote [" << tokenizedRestriction << endl;
                 exit(EXIT_FAILURE);
               }
           }
@@ -3314,25 +3304,22 @@ MarkovSwitchingStatement::writeCOutput(ostream &output, const string &basename)
   using namespace boost;
   vector<string> tokenizedDomain;
   split(tokenizedDomain, it->second, is_any_of("[ ]"), token_compress_on);
-  for (vector<string>::iterator itvs = tokenizedDomain.begin();
-       itvs != tokenizedDomain.end(); itvs++)
-    if (!itvs->empty())
-      output << "duration.push_back(" << *itvs << ");" << endl;
+  for (auto & itvs : tokenizedDomain)
+    if (!itvs.empty())
+      output << "duration.push_back(" << itvs << ");" << endl;
 
   OptionsList::symbol_list_options_t::const_iterator itsl
     = options_list.symbol_list_options.find("ms.parameters");
   assert(itsl != options_list.symbol_list_options.end());
   vector<string> parameters = itsl->second.get_symbols();
   output << "parameters.clear();" << endl;
-  for (vector<string>::iterator itp = parameters.begin();
-       itp != parameters.end(); itp++)
-    output << "parameters.push_back(param_names[\"" << *itp << "\"]);" << endl;
+  for (auto & parameter : parameters)
+    output << "parameters.push_back(param_names[\"" << parameter << "\"]);" << endl;
 
   output << "restriction_map.clear();" << endl;
-  for (map <pair<int, int >, double >::iterator itrm = restriction_map.begin();
-       itrm != restriction_map.end(); itrm++)
-    output << "restriction_map[make_pair(" << itrm->first.first << ","
-           << itrm->first.second << ")] = " << itrm->second << ";" << endl;
+  for (auto & itrm : restriction_map)
+    output << "restriction_map[make_pair(" << itrm.first.first << ","
+           << itrm.first.second << ")] = " << itrm.second << ";" << endl;
 
   output << "msdsgeinfo->addMarkovSwitching(new MarkovSwitching(" << endl
          << "     chain, number_of_regimes, number_of_lags, number_of_lags_was_passed, parameters, duration, restriction_map));" << endl;
@@ -3417,9 +3404,8 @@ SvarStatement::writeOutput(ostream &output, const string &basename, bool minimal
       if (itv->second.size() > 1)
         {
           output << "[";
-          for (vector<int>::const_iterator viit = itv->second.begin();
-               viit != itv->second.end(); viit++)
-            output << *viit << ";";
+          for (int viit : itv->second)
+            output << viit << ";";
           output << "];" << endl;
         }
       else
@@ -3751,14 +3737,14 @@ JointPriorStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsoli
 void
 JointPriorStatement::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const
 {
-  for (vector<string>::const_iterator it = joint_parameters.begin(); it != joint_parameters.end(); it++)
+  for (const auto & joint_parameter : joint_parameters)
     output << "eifind = get_new_or_existing_ei_index('joint_parameter_prior_index', '"
-           << *it << "', '');" << endl
-           << "estimation_info.joint_parameter_prior_index(eifind) = {'" << *it << "'};" << endl;
+           << joint_parameter << "', '');" << endl
+           << "estimation_info.joint_parameter_prior_index(eifind) = {'" << joint_parameter << "'};" << endl;
 
   output << "key = {[";
-  for (vector<string>::const_iterator it = joint_parameters.begin(); it != joint_parameters.end(); it++)
-    output << "get_new_or_existing_ei_index('joint_parameter_prior_index', '" << *it << "', '') ..."
+  for (const auto & joint_parameter : joint_parameters)
+    output << "get_new_or_existing_ei_index('joint_parameter_prior_index', '" << joint_parameter << "', '') ..."
            << endl << "    ";
   output << "]};" << endl;
 
@@ -4023,10 +4009,9 @@ BasicPriorStatement::writeCDomain(ostream &output) const
       using namespace boost;
       vector<string> tokenizedDomain;
       split(tokenizedDomain, it_num->second, is_any_of("[ ]"), token_compress_on);
-      for (vector<string>::iterator it = tokenizedDomain.begin();
-           it != tokenizedDomain.end(); it++)
-        if (!it->empty())
-          output << "domain.push_back(" << *it << ");" << endl;
+      for (auto & it : tokenizedDomain)
+        if (!it.empty())
+          output << "domain.push_back(" << it << ");" << endl;
     }
 }
 
@@ -4841,10 +4826,9 @@ ExtendedPathStatement::writeOutput(ostream &output, const string &basename, bool
 {
   // Beware: options do not have the same name in the interface and in the M code...
 
-  for (OptionsList::num_options_t::const_iterator it = options_list.num_options.begin();
-       it != options_list.num_options.end(); ++it)
-    if (it->first != string("periods"))
-      output << "options_." << it->first << " = " << it->second << ";" << endl;
+  for (const auto & num_option : options_list.num_options)
+    if (num_option.first != string("periods"))
+      output << "options_." << num_option.first << " = " << num_option.second << ";" << endl;
 
   output << "extended_path([], " << options_list.num_options.find("periods")->second
          << ", [], options_, M_, oo_);" << endl;
@@ -4986,9 +4970,8 @@ GenerateIRFsStatement::writeOutput(ostream &output, const string &basename, bool
     return;
 
   output << "options_.irf_opt.irf_shock_graphtitles = { ";
-  for (vector<string>::const_iterator it = generate_irf_names.begin();
-       it != generate_irf_names.end(); it++)
-    output << "'" << *it << "'; ";
+  for (const auto & generate_irf_name : generate_irf_names)
+    output << "'" << generate_irf_name << "'; ";
   output << "};" << endl;
 
   output << "options_.irf_opt.irf_shocks = zeros(M_.exo_nbr, "
diff --git a/src/ConfigFile.cc b/src/ConfigFile.cc
index 64d899e1..3a80bc43 100644
--- a/src/ConfigFile.cc
+++ b/src/ConfigFile.cc
@@ -269,12 +269,11 @@ ConfigFile::getConfigFileInfo(const string &config_file)
                 {
                   vector<string> tokenizedPath;
                   split(tokenizedPath, tokenizedLine.back(), is_any_of(":"), token_compress_on);
-                  for (vector<string>::iterator it = tokenizedPath.begin();
-                       it != tokenizedPath.end(); it++)
-                    if (!it->empty())
+                  for (auto & it : tokenizedPath)
+                    if (!it.empty())
                       {
-                        trim(*it);
-                        includepath.push_back(*it);
+                        trim(it);
+                        includepath.push_back(it);
                       }
                 }
               else
@@ -524,11 +523,11 @@ void
 ConfigFile::checkPass(WarningConsolidation &warnings) const
 {
   bool global_init_file_declared = false;
-  for (vector<Hook *>::const_iterator it = hooks.begin(); it != hooks.end(); it++)
+  for (auto hook : hooks)
     {
-      const map <string, string> hookmap = (*it)->get_hooks();
-      for (map <string, string>::const_iterator mapit = hookmap.begin(); mapit != hookmap.end(); mapit++)
-        if (mapit->first.compare("global_init_file") == 0)
+      const map <string, string> hookmap = hook->get_hooks();
+      for (const auto & mapit : hookmap)
+        if (mapit.first.compare("global_init_file") == 0)
           if (global_init_file_declared == true)
             {
               cerr << "ERROR: Only one global initialization file may be provided." << endl;
@@ -548,56 +547,55 @@ ConfigFile::checkPass(WarningConsolidation &warnings) const
       exit(EXIT_FAILURE);
     }
 
-  for (map<string, SlaveNode *>::const_iterator it = slave_nodes.begin();
-       it != slave_nodes.end(); it++)
+  for (const auto & slave_node : slave_nodes)
     {
 #if !defined(_WIN32) && !defined(__CYGWIN32__)
       //For Linux/Mac, check that cpuNbr starts at 0
-      if (it->second->minCpuNbr != 0)
+      if (slave_node.second->minCpuNbr != 0)
         warnings << "WARNING: On Unix-based operating systems, you cannot specify the CPU that is "
                  << "used in parallel processing. This will be adjusted for you such that the "
                  << "same number of CPUs are used." << endl;
 #endif
-      if (!it->second->port.empty())
+      if (!slave_node.second->port.empty())
         try
           {
-            boost::lexical_cast< int >(it->second->port);
+            boost::lexical_cast< int >(slave_node.second->port);
           }
         catch (const boost::bad_lexical_cast &)
           {
-            cerr << "ERROR (node " << it->first << "): the port must be an integer." << endl;
+            cerr << "ERROR (node " << slave_node.first << "): the port must be an integer." << endl;
             exit(EXIT_FAILURE);
           }
-      if (!it->second->computerName.compare("localhost")) // We are working locally
+      if (!slave_node.second->computerName.compare("localhost")) // We are working locally
         {
-          if (!it->second->remoteDrive.empty())
+          if (!slave_node.second->remoteDrive.empty())
             {
-              cerr << "ERROR (node " << it->first << "): the RemoteDrive option may not be passed for a local node." << endl;
+              cerr << "ERROR (node " << slave_node.first << "): the RemoteDrive option may not be passed for a local node." << endl;
               exit(EXIT_FAILURE);
             }
-          if (!it->second->remoteDirectory.empty())
+          if (!slave_node.second->remoteDirectory.empty())
             {
-              cerr << "ERROR (node " << it->first << "): the RemoteDirectory option may not be passed for a local node." << endl;
+              cerr << "ERROR (node " << slave_node.first << "): the RemoteDirectory option may not be passed for a local node." << endl;
               exit(EXIT_FAILURE);
             }
         }
       else
         {
-          if (it->second->userName.empty())
+          if (slave_node.second->userName.empty())
             {
-              cerr << "ERROR (node " << it->first << "): the UserName option must be passed for every remote node." << endl;
+              cerr << "ERROR (node " << slave_node.first << "): the UserName option must be passed for every remote node." << endl;
               exit(EXIT_FAILURE);
             }
-          if (it->second->operatingSystem.compare("windows") == 0)
+          if (slave_node.second->operatingSystem.compare("windows") == 0)
             {
-              if (it->second->password.empty())
+              if (slave_node.second->password.empty())
                 {
-                  cerr << "ERROR (node " << it->first << "): the Password option must be passed under Windows for every remote node." << endl;
+                  cerr << "ERROR (node " << slave_node.first << "): the Password option must be passed under Windows for every remote node." << endl;
                   exit(EXIT_FAILURE);
                 }
-              if (it->second->remoteDrive.empty())
+              if (slave_node.second->remoteDrive.empty())
                 {
-                  cerr << "ERROR (node " << it->first << "): the RemoteDrive option must be passed under Windows for every remote node." << endl;
+                  cerr << "ERROR (node " << slave_node.first << "): the RemoteDrive option must be passed under Windows for every remote node." << endl;
                   exit(EXIT_FAILURE);
                 }
             }
@@ -616,9 +614,9 @@ ConfigFile::checkPass(WarningConsolidation &warnings) const
                 }
             }
 #endif
-          if (it->second->remoteDirectory.empty())
+          if (slave_node.second->remoteDirectory.empty())
             {
-              cerr << "ERROR (node " << it->first << "): the RemoteDirectory must be specified for every remote node." << endl;
+              cerr << "ERROR (node " << slave_node.first << "): the RemoteDirectory must be specified for every remote node." << endl;
               exit(EXIT_FAILURE);
             }
         }
@@ -637,13 +635,12 @@ ConfigFile::checkPass(WarningConsolidation &warnings) const
       exit(EXIT_FAILURE);
     }
 
-  for (map<string, Cluster *>::const_iterator it = clusters.begin();
-       it != clusters.end(); it++)
-    for (member_nodes_t::const_iterator itmn = it->second->member_nodes.begin();
-         itmn != it->second->member_nodes.end(); itmn++)
+  for (const auto & cluster : clusters)
+    for (member_nodes_t::const_iterator itmn = cluster.second->member_nodes.begin();
+         itmn != cluster.second->member_nodes.end(); itmn++)
       if (slave_nodes.find(itmn->first) == slave_nodes.end())
         {
-          cerr << "Error: node " << itmn->first << " specified in cluster " << it->first << " was not found" << endl;
+          cerr << "Error: node " << itmn->first << " specified in cluster " << cluster.first << " was not found" << endl;
           exit(EXIT_FAILURE);
         }
 }
@@ -676,21 +673,20 @@ ConfigFile::transformPass()
        it != cluster_it->second->member_nodes.end(); it++)
     weight_denominator += it->second;
 
-  for (member_nodes_t::iterator it = cluster_it->second->member_nodes.begin();
-       it != cluster_it->second->member_nodes.end(); it++)
-    it->second /= weight_denominator;
+  for (auto & member_node : cluster_it->second->member_nodes)
+    member_node.second /= weight_denominator;
 }
 
 vector<string>
 ConfigFile::getIncludePaths() const
 {
   vector<string> include_paths;
-  for (vector<Path *>::const_iterator it = paths.begin(); it != paths.end(); it++)
+  for (auto path : paths)
     {
-      map <string, vector<string> > pathmap = (*it)->get_paths();
+      map <string, vector<string> > pathmap = path->get_paths();
       for (map <string, vector<string> >::const_iterator mapit = pathmap.begin(); mapit != pathmap.end(); mapit++)
-        for (vector<string>::const_iterator vecit = mapit->second.begin(); vecit != mapit->second.end(); vecit++)
-          include_paths.push_back(*vecit);
+        for (const auto & vecit : mapit->second)
+          include_paths.push_back(vecit);
     }
   return include_paths;
 }
@@ -698,9 +694,9 @@ ConfigFile::getIncludePaths() const
 void
 ConfigFile::writeHooks(ostream &output) const
 {
-  for (vector<Hook *>::const_iterator it = hooks.begin(); it != hooks.end(); it++)
+  for (auto hook : hooks)
     {
-      map <string, string> hookmap = (*it)->get_hooks();
+      map <string, string> hookmap = hook->get_hooks();
       for (map <string, string>::const_iterator mapit = hookmap.begin(); mapit != hookmap.end(); mapit++)
         output << "options_." << mapit->first << " = '" << mapit->second << "';" << endl;
     }
@@ -719,13 +715,12 @@ ConfigFile::writeCluster(ostream &output) const
     cluster_it = clusters.find(cluster_name);
 
   int i = 1;
-  for (map<string, SlaveNode *>::const_iterator it = slave_nodes.begin();
-       it != slave_nodes.end(); it++)
+  for (const auto & slave_node : slave_nodes)
     {
       bool slave_node_in_member_nodes = false;
       for (member_nodes_t::const_iterator itmn = cluster_it->second->member_nodes.begin();
            itmn != cluster_it->second->member_nodes.end(); itmn++)
-        if (!it->first.compare(itmn->first))
+        if (!slave_node.first.compare(itmn->first))
           slave_node_in_member_nodes = true;
 
       if (!slave_node_in_member_nodes)
@@ -736,25 +731,25 @@ ConfigFile::writeCluster(ostream &output) const
         output << "(" << i << ")";
       i++;
       output << " = struct('Local', ";
-      if (it->second->computerName.compare("localhost"))
+      if (slave_node.second->computerName.compare("localhost"))
         output << "0, ";
       else
         output << "1, ";
 
-      output << "'ComputerName', '" << it->second->computerName << "', "
-             << "'Port', '" << it->second->port << "', "
-             << "'CPUnbr', [" << it->second->minCpuNbr << ":" << it->second->maxCpuNbr << "], "
-             << "'UserName', '" << it->second->userName << "', "
-             << "'Password', '" << it->second->password << "', "
-             << "'RemoteDrive', '" << it->second->remoteDrive << "', "
-             << "'RemoteDirectory', '" << it->second->remoteDirectory << "', "
-             << "'DynarePath', '" << it->second->dynarePath << "', "
-             << "'MatlabOctavePath', '" << it->second->matlabOctavePath << "', "
-             << "'OperatingSystem', '" << it->second->operatingSystem << "', "
-             << "'NodeWeight', '" << (cluster_it->second->member_nodes.find(it->first))->second << "', "
-             << "'NumberOfThreadsPerJob', " << it->second->numberOfThreadsPerJob << ", ";
+      output << "'ComputerName', '" << slave_node.second->computerName << "', "
+             << "'Port', '" << slave_node.second->port << "', "
+             << "'CPUnbr', [" << slave_node.second->minCpuNbr << ":" << slave_node.second->maxCpuNbr << "], "
+             << "'UserName', '" << slave_node.second->userName << "', "
+             << "'Password', '" << slave_node.second->password << "', "
+             << "'RemoteDrive', '" << slave_node.second->remoteDrive << "', "
+             << "'RemoteDirectory', '" << slave_node.second->remoteDirectory << "', "
+             << "'DynarePath', '" << slave_node.second->dynarePath << "', "
+             << "'MatlabOctavePath', '" << slave_node.second->matlabOctavePath << "', "
+             << "'OperatingSystem', '" << slave_node.second->operatingSystem << "', "
+             << "'NodeWeight', '" << (cluster_it->second->member_nodes.find(slave_node.first))->second << "', "
+             << "'NumberOfThreadsPerJob', " << slave_node.second->numberOfThreadsPerJob << ", ";
 
-      if (it->second->singleCompThread)
+      if (slave_node.second->singleCompThread)
         output << "'SingleCompThread', 'true');" << endl;
       else
         output << "'SingleCompThread', 'false');" << endl;
diff --git a/src/DataTree.cc b/src/DataTree.cc
index f193d0da..2c5c9c79 100644
--- a/src/DataTree.cc
+++ b/src/DataTree.cc
@@ -46,8 +46,8 @@ DataTree::DataTree(SymbolTable &symbol_table_arg,
 
 DataTree::~DataTree()
 {
-  for (node_list_t::iterator it = node_list.begin(); it != node_list.end(); it++)
-    delete *it;
+  for (auto & it : node_list)
+    delete it;
 }
 
 expr_t
@@ -75,9 +75,8 @@ DataTree::AddVariableInternal(int symb_id, int lag)
 bool
 DataTree::ParamUsedWithLeadLagInternal() const
 {
-  for (variable_node_map_t::const_iterator it = variable_node_map.begin();
-       it != variable_node_map.end(); it++)
-    if (symbol_table.getType(it->first.first) == eParameter && it->first.second != 0)
+  for (const auto & it : variable_node_map)
+    if (symbol_table.getType(it.first.first) == eParameter && it.first.second != 0)
       return true;
   return false;
 }
@@ -588,9 +587,8 @@ DataTree::AddSecondDerivExternalFunction(int top_level_symb_id, const vector<exp
 bool
 DataTree::isSymbolUsed(int symb_id) const
 {
-  for (variable_node_map_t::const_iterator it = variable_node_map.begin();
-       it != variable_node_map.end(); it++)
-    if (it->first.first == symb_id)
+  for (const auto & it : variable_node_map)
+    if (it.first.first == symb_id)
       return true;
 
   if (local_variables_table.find(symb_id) != local_variables_table.end())
@@ -637,9 +635,8 @@ DataTree::getDynJacobianCol(int deriv_id) const throw (UnknownDerivIDException)
 bool
 DataTree::isUnaryOpUsed(UnaryOpcode opcode) const
 {
-  for (unary_op_node_map_t::const_iterator it = unary_op_node_map.begin();
-       it != unary_op_node_map.end(); it++)
-    if (it->first.first.second == opcode)
+  for (const auto & it : unary_op_node_map)
+    if (it.first.first.second == opcode)
       return true;
 
   return false;
@@ -648,9 +645,8 @@ DataTree::isUnaryOpUsed(UnaryOpcode opcode) const
 bool
 DataTree::isBinaryOpUsed(BinaryOpcode opcode) const
 {
-  for (binary_op_node_map_t::const_iterator it = binary_op_node_map.begin();
-       it != binary_op_node_map.end(); it++)
-    if (it->first.second == opcode)
+  for (const auto & it : binary_op_node_map)
+    if (it.first.second == opcode)
       return true;
 
   return false;
@@ -659,9 +655,8 @@ DataTree::isBinaryOpUsed(BinaryOpcode opcode) const
 bool
 DataTree::isTrinaryOpUsed(TrinaryOpcode opcode) const
 {
-  for (trinary_op_node_map_t::const_iterator it = trinary_op_node_map.begin();
-       it != trinary_op_node_map.end(); it++)
-    if (it->first.second == opcode)
+  for (const auto & it : trinary_op_node_map)
+    if (it.first.second == opcode)
       return true;
 
   return false;
@@ -670,9 +665,8 @@ DataTree::isTrinaryOpUsed(TrinaryOpcode opcode) const
 bool
 DataTree::isExternalFunctionUsed(int symb_id) const
 {
-  for (external_function_node_map_t::const_iterator it = external_function_node_map.begin();
-       it != external_function_node_map.end(); it++)
-    if (it->first.second == symb_id)
+  for (const auto & it : external_function_node_map)
+    if (it.first.second == symb_id)
       return true;
 
   return false;
@@ -681,9 +675,8 @@ DataTree::isExternalFunctionUsed(int symb_id) const
 bool
 DataTree::isFirstDerivExternalFunctionUsed(int symb_id) const
 {
-  for (first_deriv_external_function_node_map_t::const_iterator it = first_deriv_external_function_node_map.begin();
-       it != first_deriv_external_function_node_map.end(); it++)
-    if (it->first.second == symb_id)
+  for (const auto & it : first_deriv_external_function_node_map)
+    if (it.first.second == symb_id)
       return true;
 
   return false;
@@ -692,9 +685,8 @@ DataTree::isFirstDerivExternalFunctionUsed(int symb_id) const
 bool
 DataTree::isSecondDerivExternalFunctionUsed(int symb_id) const
 {
-  for (second_deriv_external_function_node_map_t::const_iterator it = second_deriv_external_function_node_map.begin();
-       it != second_deriv_external_function_node_map.end(); it++)
-    if (it->first.second == symb_id)
+  for (const auto & it : second_deriv_external_function_node_map)
+    if (it.first.second == symb_id)
       return true;
 
   return false;
@@ -704,10 +696,9 @@ int
 DataTree::minLagForSymbol(int symb_id) const
 {
   int r = 0;
-  for (variable_node_map_t::const_iterator it = variable_node_map.begin();
-       it != variable_node_map.end(); ++it)
-    if (it->first.first == symb_id && it->first.second < r)
-      r = it->first.second;
+  for (const auto & it : variable_node_map)
+    if (it.first.first == symb_id && it.first.second < r)
+      r = it.first.second;
   return r;
 }
 
diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 6f1e8614..3ee782dc 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -206,9 +206,8 @@ DynamicModel::computeTemporaryTermsMapping()
 {
   // Add a mapping form node ID to temporary terms order
   int j = 0;
-  for (temporary_terms_t::const_iterator it = temporary_terms.begin();
-       it != temporary_terms.end(); it++)
-    map_idx[(*it)->idx] = j++;
+  for (auto temporary_term : temporary_terms)
+    map_idx[temporary_term->idx] = j++;
 }
 
 void
@@ -400,9 +399,8 @@ DynamicModel::writeModelEquationsOrdered_M(const string &dynamic_basename) const
       if (v_temporary_terms_inuse[block].size())
         {
           tmp_output.str("");
-          for (temporary_terms_inuse_t::const_iterator it = v_temporary_terms_inuse[block].begin();
-               it != v_temporary_terms_inuse[block].end(); it++)
-            tmp_output << " T" << *it;
+          for (int it : v_temporary_terms_inuse[block])
+            tmp_output << " T" << it;
           output << "  global" << tmp_output.str() << ";\n";
         }
       if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
@@ -415,12 +413,11 @@ DynamicModel::writeModelEquationsOrdered_M(const string &dynamic_basename) const
                 {
                   output << "  " << "% //Temporary variables initialization" << endl
                          << "  " << "T_zeros = zeros(y_kmin+periods, 1);" << endl;
-                  for (temporary_terms_t::const_iterator it = v_temporary_terms[block][i].begin();
-                       it != v_temporary_terms[block][i].end(); it++)
+                  for (auto it : v_temporary_terms[block][i])
                     {
                       output << "  ";
                       // In the following, "Static" is used to avoid getting the "(it_)" subscripting
-                      (*it)->writeOutput(output, oMatlabStaticModelSparse, local_temporary_terms, {});
+                      it->writeOutput(output, oMatlabStaticModelSparse, local_temporary_terms, {});
                       output << " = T_zeros;" << endl;
                     }
                 }
@@ -458,18 +455,17 @@ DynamicModel::writeModelEquationsOrdered_M(const string &dynamic_basename) const
           if (v_temporary_terms[block].size())
             {
               output << "  " << "% //Temporary variables" << endl;
-              for (temporary_terms_t::const_iterator it = v_temporary_terms[block][i].begin();
-                   it != v_temporary_terms[block][i].end(); it++)
+              for (auto it : v_temporary_terms[block][i])
                 {
-                  if (dynamic_cast<AbstractExternalFunctionNode *>(*it) != NULL)
-                    (*it)->writeExternalFunctionOutput(output, local_output_type, tt2, temporary_terms_idxs, tef_terms);
+                  if (dynamic_cast<AbstractExternalFunctionNode *>(it) != NULL)
+                    it->writeExternalFunctionOutput(output, local_output_type, tt2, temporary_terms_idxs, tef_terms);
 
                   output << "  " <<  sps;
-                  (*it)->writeOutput(output, local_output_type, local_temporary_terms, {}, tef_terms);
+                  it->writeOutput(output, local_output_type, local_temporary_terms, {}, tef_terms);
                   output << " = ";
-                  (*it)->writeOutput(output, local_output_type, tt2, {}, tef_terms);
+                  it->writeOutput(output, local_output_type, tt2, {}, tef_terms);
                   // Insert current node into tt2
-                  tt2.insert(*it);
+                  tt2.insert(it);
                   output << ";" << endl;
                 }
             }
@@ -836,18 +832,17 @@ DynamicModel::writeModelEquationsCode(string &file_name, const string &bin_basen
 
   map<pair< int, pair<int, int> >, expr_t> first_derivatives_reordered_endo;
   map<pair< pair<int, int>, pair<int, int> >, expr_t>  first_derivatives_reordered_exo;
-  for (first_derivatives_t::const_iterator it = first_derivatives.begin();
-       it != first_derivatives.end(); it++)
+  for (const auto & first_derivative : first_derivatives)
     {
-      int deriv_id = it->first.second;
-      unsigned int eq = it->first.first;
+      int deriv_id = first_derivative.first.second;
+      unsigned int eq = first_derivative.first.first;
       int symb = getSymbIDByDerivID(deriv_id);
       unsigned int var = symbol_table.getTypeSpecificID(symb);
       int lag = getLagByDerivID(deriv_id);
       if (getTypeByDerivID(deriv_id) == eEndogenous)
-        first_derivatives_reordered_endo[make_pair(lag, make_pair(var, eq))] = it->second;
+        first_derivatives_reordered_endo[make_pair(lag, make_pair(var, eq))] = first_derivative.second;
       else if (getTypeByDerivID(deriv_id) == eExogenous || getTypeByDerivID(deriv_id) == eExogenousDet)
-        first_derivatives_reordered_exo[make_pair(make_pair(lag, getTypeByDerivID(deriv_id)), make_pair(var, eq))] = it->second;
+        first_derivatives_reordered_exo[make_pair(make_pair(lag, getTypeByDerivID(deriv_id)), make_pair(var, eq))] = first_derivative.second;
     }
   int prev_var = -1;
   int prev_lag = -999999999;
@@ -928,14 +923,13 @@ DynamicModel::writeModelEquationsCode(string &file_name, const string &bin_basen
   vector<vector<pair<pair<int, int>, int > > > derivatives;
   derivatives.resize(symbol_table.endo_nbr());
   count_u = symbol_table.endo_nbr();
-  for (first_derivatives_t::const_iterator it = first_derivatives.begin();
-       it != first_derivatives.end(); it++)
+  for (const auto & first_derivative : first_derivatives)
     {
-      int deriv_id = it->first.second;
+      int deriv_id = first_derivative.first.second;
       if (getTypeByDerivID(deriv_id) == eEndogenous)
         {
-          expr_t d1 = it->second;
-          unsigned int eq = it->first.first;
+          expr_t d1 = first_derivative.second;
+          unsigned int eq = first_derivative.first.first;
           int symb = getSymbIDByDerivID(deriv_id);
           unsigned int var = symbol_table.getTypeSpecificID(symb);
           int lag = getLagByDerivID(deriv_id);
@@ -1142,8 +1136,8 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
         }
       unsigned int count_col_det_exo = 0;
       vector<unsigned int> exo_det;
-      for (lag_var_t::const_iterator it = exo_det_block[block].begin(); it != exo_det_block[block].end(); it++)
-        for (var_t::const_iterator it1 = it->second.begin(); it1 != it->second.end(); it1++)
+      for (const auto & it : exo_det_block[block])
+        for (var_t::const_iterator it1 = it.second.begin(); it1 != it.second.end(); it1++)
           {
             count_col_det_exo++;
             if (find(exo_det.begin(), exo_det.end(), *it1) == exo_det.end())
@@ -1152,8 +1146,8 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
 
       unsigned int count_col_exo = 0;
       vector<unsigned int> exo;
-      for (lag_var_t::const_iterator it = exo_block[block].begin(); it != exo_block[block].end(); it++)
-        for (var_t::const_iterator it1 = it->second.begin(); it1 != it->second.end(); it1++)
+      for (const auto & it : exo_block[block])
+        for (var_t::const_iterator it1 = it.second.begin(); it1 != it.second.end(); it1++)
           {
             count_col_exo++;
             if (find(exo.begin(), exo.end(), *it1) == exo.end())
@@ -1162,8 +1156,8 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
 
       vector<unsigned int> other_endo;
       unsigned int count_col_other_endo = 0;
-      for (lag_var_t::const_iterator it = other_endo_block[block].begin(); it != other_endo_block[block].end(); it++)
-        for (var_t::const_iterator it1 = it->second.begin(); it1 != it->second.end(); it1++)
+      for (const auto & it : other_endo_block[block])
+        for (var_t::const_iterator it1 = it.second.begin(); it1 != it.second.end(); it1++)
           {
             count_col_other_endo++;
             if (find(other_endo.begin(), other_endo.end(), *it1) == other_endo.end())
@@ -1202,19 +1196,18 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
           tt2.clear();
           if (v_temporary_terms[block][i].size())
             {
-              for (temporary_terms_t::const_iterator it = v_temporary_terms[block][i].begin();
-                   it != v_temporary_terms[block][i].end(); it++)
+              for (auto it : v_temporary_terms[block][i])
                 {
-                  if (dynamic_cast<AbstractExternalFunctionNode *>(*it) != NULL)
-                    (*it)->compileExternalFunctionOutput(code_file, instruction_number, false, tt2, map_idx, true, false, tef_terms);
+                  if (dynamic_cast<AbstractExternalFunctionNode *>(it) != NULL)
+                    it->compileExternalFunctionOutput(code_file, instruction_number, false, tt2, map_idx, true, false, tef_terms);
 
-                  FNUMEXPR_ fnumexpr(TemporaryTerm, (int)(map_idx.find((*it)->idx)->second));
+                  FNUMEXPR_ fnumexpr(TemporaryTerm, (int)(map_idx.find(it->idx)->second));
                   fnumexpr.write(code_file, instruction_number);
-                  (*it)->compile(code_file, instruction_number, false, tt2, map_idx, true, false, tef_terms);
-                  FSTPT_ fstpt((int)(map_idx.find((*it)->idx)->second));
+                  it->compile(code_file, instruction_number, false, tt2, map_idx, true, false, tef_terms);
+                  FSTPT_ fstpt((int)(map_idx.find(it->idx)->second));
                   fstpt.write(code_file, instruction_number);
                   // Insert current node into tt2
-                  tt2.insert(*it);
+                  tt2.insert(it);
 #ifdef DEBUGC
                   cout << "FSTPT " << v << "\n";
                   instruction_number++;
@@ -1544,16 +1537,15 @@ DynamicModel::getVarExpectationFunctionsToWrite() const
 void
 DynamicModel::writeVarExpectationCalls(ostream &output) const
 {
-  for (map<string, set<int> >::const_iterator it = var_expectation_functions_to_write.begin();
-       it != var_expectation_functions_to_write.end(); it++)
+  for (const auto & it : var_expectation_functions_to_write)
     {
       int i = 0;
-      output << "dynamic_var_forecast_" << it->first << " = "
-             << "var_forecast_" << it->first << "(y);" << endl;
+      output << "dynamic_var_forecast_" << it.first << " = "
+             << "var_forecast_" << it.first << "(y);" << endl;
 
-      for (set<int>::const_iterator it1 = it->second.begin(); it1 != it->second.end(); it1++)
-        output << "dynamic_var_forecast_" << it->first << "_" << *it1 << " = "
-               << "dynamic_var_forecast_" << it->first << "(" << ++i << ", :);" << endl;
+      for (set<int>::const_iterator it1 = it.second.begin(); it1 != it.second.end(); it1++)
+        output << "dynamic_var_forecast_" << it.first << "_" << *it1 << " = "
+               << "dynamic_var_forecast_" << it.first << "(" << ++i << ", :);" << endl;
     }
 }
 
@@ -1737,11 +1729,11 @@ DynamicModel::setNonZeroHessianEquations(map<int, string> &eqs)
     if (nonzero_hessian_eqs.find(it->first.first) == nonzero_hessian_eqs.end())
       {
         nonzero_hessian_eqs[it->first.first] = "";
-        for (size_t i = 0; i < equation_tags.size(); i++)
-          if (equation_tags[i].first == it->first.first)
-            if (equation_tags[i].second.first == "name")
+        for (auto & equation_tag : equation_tags)
+          if (equation_tag.first == it->first.first)
+            if (equation_tag.second.first == "name")
               {
-                nonzero_hessian_eqs[it->first.first] = equation_tags[i].second.second;
+                nonzero_hessian_eqs[it->first.first] = equation_tag.second.second;
                 break;
               }
       }
@@ -1832,27 +1824,25 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri
   //Temporary variables declaration
   OK = true;
   ostringstream tmp_output;
-  for (temporary_terms_t::const_iterator it = temporary_terms.begin();
-       it != temporary_terms.end(); it++)
+  for (auto temporary_term : temporary_terms)
     {
       if (OK)
         OK = false;
       else
         tmp_output << " ";
       // In the following, "Static" is used to avoid getting the "(it_)" subscripting
-      (*it)->writeOutput(tmp_output, oMatlabStaticModelSparse, temporary_terms, {});
+      temporary_term->writeOutput(tmp_output, oMatlabStaticModelSparse, temporary_terms, {});
     }
   if (tmp_output.str().length() > 0)
     mDynamicModelFile << "  global " << tmp_output.str() << ";\n";
 
   mDynamicModelFile << "  T_init=zeros(1,options_.periods+M_.maximum_lag+M_.maximum_lead);\n";
   tmp_output.str("");
-  for (temporary_terms_t::const_iterator it = temporary_terms.begin();
-       it != temporary_terms.end(); it++)
+  for (auto temporary_term : temporary_terms)
     {
       tmp_output << "  ";
       // In the following, "Static" is used to avoid getting the "(it_)" subscripting
-      (*it)->writeOutput(tmp_output, oMatlabStaticModelSparse, temporary_terms, {});
+      temporary_term->writeOutput(tmp_output, oMatlabStaticModelSparse, temporary_terms, {});
       tmp_output << "=T_init;\n";
     }
   if (tmp_output.str().length() > 0)
@@ -2343,9 +2333,8 @@ DynamicModel::writeDynamicModel(const string &dynamic_basename, ostream &Dynamic
   deriv_node_temp_terms_t tef_terms;
   temporary_terms_t temp_term_union;
 
-  for (map<expr_t, expr_t>::const_iterator it = temporary_terms_mlv.begin();
-       it != temporary_terms_mlv.end(); it++)
-    temp_term_union.insert(it->first);
+  for (auto it : temporary_terms_mlv)
+    temp_term_union.insert(it.first);
   writeModelLocalVariableTemporaryTerms(temp_term_union, temporary_terms_mlv,
                                         model_tt_output, output_type, tef_terms);
 
@@ -2369,12 +2358,11 @@ DynamicModel::writeDynamicModel(const string &dynamic_basename, ostream &Dynamic
                           jacobian_tt_output, output_type, tef_terms);
       temp_term_union.insert(temporary_terms_g1.begin(), temporary_terms_g1.end());
 
-      for (first_derivatives_t::const_iterator it = first_derivatives.begin();
-           it != first_derivatives.end(); it++)
+      for (const auto & first_derivative : first_derivatives)
         {
-          int eq = it->first.first;
-          int var = it->first.second;
-          expr_t d1 = it->second;
+          int eq = first_derivative.first.first;
+          int var = first_derivative.first.second;
+          expr_t d1 = first_derivative.second;
 
           jacobianHelper(jacobian_output, eq, getDynJacobianCol(var), output_type);
           jacobian_output << "=";
@@ -2394,13 +2382,12 @@ DynamicModel::writeDynamicModel(const string &dynamic_basename, ostream &Dynamic
       temp_term_union.insert(temporary_terms_g2.begin(), temporary_terms_g2.end());
 
       int k = 0; // Keep the line of a 2nd derivative in v2
-      for (second_derivatives_t::const_iterator it = second_derivatives.begin();
-           it != second_derivatives.end(); it++)
+      for (const auto & second_derivative : second_derivatives)
         {
-          int eq = it->first.first;
-          int var1 = it->first.second.first;
-          int var2 = it->first.second.second;
-          expr_t d2 = it->second;
+          int eq = second_derivative.first.first;
+          int var1 = second_derivative.first.second.first;
+          int var2 = second_derivative.first.second.second;
+          expr_t d2 = second_derivative.second;
 
           int id1 = getDynJacobianCol(var1);
           int id2 = getDynJacobianCol(var2);
@@ -2465,14 +2452,13 @@ DynamicModel::writeDynamicModel(const string &dynamic_basename, ostream &Dynamic
       temp_term_union.insert(temporary_terms_g3.begin(), temporary_terms_g3.end());
 
       int k = 0; // Keep the line of a 3rd derivative in v3
-      for (third_derivatives_t::const_iterator it = third_derivatives.begin();
-           it != third_derivatives.end(); it++)
+      for (const auto & third_derivative : third_derivatives)
         {
-          int eq = it->first.first;
-          int var1 = it->first.second.first;
-          int var2 = it->first.second.second.first;
-          int var3 = it->first.second.second.second;
-          expr_t d3 = it->second;
+          int eq = third_derivative.first.first;
+          int var1 = third_derivative.first.second.first;
+          int var2 = third_derivative.first.second.second.first;
+          int var3 = third_derivative.first.second.second.second;
+          expr_t d3 = third_derivative.second;
 
           int id1 = getDynJacobianCol(var1);
           int id2 = getDynJacobianCol(var2);
@@ -2513,10 +2499,10 @@ DynamicModel::writeDynamicModel(const string &dynamic_basename, ostream &Dynamic
           cols.insert(id3 * hessianColsNbr + id2 * dynJacobianColsNbr + id1);
 
           int k2 = 1; // Keeps the offset of the permutation relative to k
-          for (set<int>::iterator it2 = cols.begin(); it2 != cols.end(); it2++)
-            if (*it2 != ref_col)
+          for (int col : cols)
+            if (col != ref_col)
               if (output_type == oJuliaDynamicModel)
-                third_derivatives_output << "    @inbounds g3[" << eq + 1 << "," << *it2 + 1 << "] = "
+                third_derivatives_output << "    @inbounds g3[" << eq + 1 << "," << col + 1 << "] = "
                                          << for_sym.str() << endl;
               else
                 {
@@ -2524,7 +2510,7 @@ DynamicModel::writeDynamicModel(const string &dynamic_basename, ostream &Dynamic
                   third_derivatives_output << "=" << eq + 1 << ";" << endl;
 
                   sparseHelper(3, third_derivatives_output, k+k2, 1, output_type);
-                  third_derivatives_output << "=" << *it2 + 1 << ";" << endl;
+                  third_derivatives_output << "=" << col + 1 << ";" << endl;
 
                   sparseHelper(3, third_derivatives_output, k+k2, 2, output_type);
                   third_derivatives_output << "=";
@@ -2961,20 +2947,20 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
   if (julia)
     {
       output << modstruct << "equation_tags = [" << endl;
-      for (size_t i = 0; i < equation_tags.size(); i++)
+      for (const auto & equation_tag : equation_tags)
         output << "                       EquationTag("
-               << equation_tags[i].first + 1 << " , \""
-               << equation_tags[i].second.first << "\" , \""
-               << equation_tags[i].second.second << "\")" << endl;
+               << equation_tag.first + 1 << " , \""
+               << equation_tag.second.first << "\" , \""
+               << equation_tag.second.second << "\")" << endl;
       output << "                      ]" << endl;
     }
   else
     {
       output << modstruct << "equations_tags = {" << endl;
-      for (size_t i = 0; i < equation_tags.size(); i++)
-        output << "  " << equation_tags[i].first + 1 << " , '"
-               << equation_tags[i].second.first << "' , '"
-               << equation_tags[i].second.second << "' ;" << endl;
+      for (const auto & equation_tag : equation_tags)
+        output << "  " << equation_tag.first + 1 << " , '"
+               << equation_tag.second.first << "' , '"
+               << equation_tag.second.second << "' ;" << endl;
       output << "};" << endl;
     }
 
@@ -3031,18 +3017,18 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
             }
           set<int> exogenous;
           exogenous.clear();
-          for (lag_var_t::const_iterator it = exo_block[block].begin(); it != exo_block[block].end(); it++)
-            for (var_t::const_iterator it1 = it->second.begin(); it1 != it->second.end(); it1++)
+          for (const auto & it : exo_block[block])
+            for (var_t::const_iterator it1 = it.second.begin(); it1 != it.second.end(); it1++)
               exogenous.insert(*it1);
           set<int> exogenous_det;
           exogenous_det.clear();
-          for (lag_var_t::const_iterator it = exo_det_block[block].begin(); it != exo_det_block[block].end(); it++)
-            for (var_t::const_iterator it1 = it->second.begin(); it1 != it->second.end(); it1++)
+          for (const auto & it : exo_det_block[block])
+            for (var_t::const_iterator it1 = it.second.begin(); it1 != it.second.end(); it1++)
               exogenous_det.insert(*it1);
           set<int> other_endogenous;
           other_endogenous.clear();
-          for (lag_var_t::const_iterator it = other_endo_block[block].begin(); it != other_endo_block[block].end(); it++)
-            for (var_t::const_iterator it1 = it->second.begin(); it1 != it->second.end(); it1++)
+          for (const auto & it : other_endo_block[block])
+            for (var_t::const_iterator it1 = it.second.begin(); it1 != it.second.end(); it1++)
               other_endogenous.insert(*it1);
           output << "block_structure.block(" << block+1 << ").Simulation_Type = " << simulation_type << ";\n";
           output << "block_structure.block(" << block+1 << ").maximum_lag = " << max_lag << ";\n";
@@ -3060,20 +3046,20 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
           output << "block_structure.block(" << block+1 << ").exo_nbr = " << getBlockExoSize(block) << ";\n";
           output << "block_structure.block(" << block+1 << ").exogenous = [";
           int i = 0;
-          for (set<int>::iterator it_exogenous = exogenous.begin(); it_exogenous != exogenous.end(); it_exogenous++)
-            if (*it_exogenous >= 0)
+          for (int exogenou : exogenous)
+            if (exogenou >= 0)
               {
-                output << " " << *it_exogenous+1;
+                output << " " << exogenou+1;
                 i++;
               }
           output << "];\n";
 
           output << "block_structure.block(" << block+1 << ").exogenous_det = [";
           i = 0;
-          for (set<int>::iterator it_exogenous_det = exogenous_det.begin(); it_exogenous_det != exogenous_det.end(); it_exogenous_det++)
-            if (*it_exogenous_det >= 0)
+          for (int it_exogenous_det : exogenous_det)
+            if (it_exogenous_det >= 0)
               {
-                output << " " << *it_exogenous_det+1;
+                output << " " << it_exogenous_det+1;
                 i++;
               }
           output << "];\n";
@@ -3081,17 +3067,17 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
 
           output << "block_structure.block(" << block+1 << ").other_endogenous = [";
           i = 0;
-          for (set<int>::iterator it_other_endogenous = other_endogenous.begin(); it_other_endogenous != other_endogenous.end(); it_other_endogenous++)
-            if (*it_other_endogenous >= 0)
+          for (int other_endogenou : other_endogenous)
+            if (other_endogenou >= 0)
               {
-                output << " " << *it_other_endogenous+1;
+                output << " " << other_endogenou+1;
                 i++;
               }
           output << "];\n";
           output << "block_structure.block(" << block+1 << ").other_endogenous_block = [";
           i = 0;
-          for (set<int>::iterator it_other_endogenous = other_endogenous.begin(); it_other_endogenous != other_endogenous.end(); it_other_endogenous++)
-            if (*it_other_endogenous >= 0)
+          for (int other_endogenou : other_endogenous)
+            if (other_endogenou >= 0)
               {
                 bool OK = true;
                 unsigned int j;
@@ -3099,7 +3085,7 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
                   for (unsigned int k = 0; k < getBlockSize(j) && OK; k++)
                     {
                       //printf("*it_other_endogenous=%d, getBlockVariableID(%d, %d)=%d\n",*it_other_endogenous, j, k, getBlockVariableID(j, k));
-                      OK = *it_other_endogenous != getBlockVariableID(j, k);
+                      OK = other_endogenou != getBlockVariableID(j, k);
                     }
                 if (!OK)
                   output << " " << j;
@@ -3110,12 +3096,12 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
           //vector<int> inter_state_var;
           output << "block_structure.block(" << block+1 << ").tm1 = zeros(" << i << ", " << state_var.size() << ");\n";
           int count_other_endogenous = 1;
-          for (set<int>::const_iterator it_other_endogenous = other_endogenous.begin(); it_other_endogenous != other_endogenous.end(); it_other_endogenous++)
+          for (int other_endogenou : other_endogenous)
             {
               for (vector<int>::const_iterator it = state_var.begin(); it != state_var.end(); it++)
                 {
                   //cout << "block = " << block+1 << " state_var = " << *it << " it_other_endogenous=" << *it_other_endogenous + 1 << "\n";
-                  if (*it == *it_other_endogenous + 1)
+                  if (*it == other_endogenou + 1)
                     {
                       output << "block_structure.block(" << block+1 << ").tm1("
                              << count_other_endogenous << ", "
@@ -3134,8 +3120,8 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
           tmp_s.str("");
           count_lead_lag_incidence = 0;
           dynamic_jacob_map_t reordered_dynamic_jacobian;
-          for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != blocks_derivatives[block].end(); it++)
-            reordered_dynamic_jacobian[make_pair(it->second.first, make_pair(it->first.second, it->first.first))] = it->second.second;
+          for (const auto & it : blocks_derivatives[block])
+            reordered_dynamic_jacobian[make_pair(it.second.first, make_pair(it.first.second, it.first.first))] = it.second.second;
           output << "block_structure.block(" << block+1 << ").lead_lag_incidence = [];\n";
           int last_var = -1;
           vector<int> local_state_var;
@@ -3203,13 +3189,13 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
           for (int lag = -1; lag <= 1; lag++)
             {
               tmp_s.str("");
-              for (set<int>::iterator it_other_endogenous = other_endogenous.begin(); it_other_endogenous != other_endogenous.end(); it_other_endogenous++)
+              for (int other_endogenou : other_endogenous)
                 {
                   bool done = false;
                   for (int i = 0; i < block_size; i++)
                     {
                       unsigned int eq = getBlockEquationID(block, i);
-                      derivative_t::const_iterator it = derivative_other_endo[block].find(make_pair(lag, make_pair(eq, *it_other_endogenous)));
+                      derivative_t::const_iterator it = derivative_other_endo[block].find(make_pair(lag, make_pair(eq, other_endogenou)));
                       if (it != derivative_other_endo[block].end())
                         {
                           count_lead_lag_incidence++;
@@ -3248,13 +3234,12 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
         state_equ.push_back(equation_reordered[variable_inv_reordered[*it - 1]]+1);
 
       map<pair< int, pair<int, int> >,  int>  lag_row_incidence;
-      for (first_derivatives_t::const_iterator it = first_derivatives.begin();
-           it != first_derivatives.end(); it++)
+      for (const auto & first_derivative : first_derivatives)
         {
-          int deriv_id = it->first.second;
+          int deriv_id = first_derivative.first.second;
           if (getTypeByDerivID(deriv_id) == eEndogenous)
             {
-              int eq = it->first.first;
+              int eq = first_derivative.first.first;
               int symb = getSymbIDByDerivID(deriv_id);
               int var = symbol_table.getTypeSpecificID(symb);
               int lag = getLagByDerivID(deriv_id);
@@ -3350,8 +3335,8 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
 
                     }
                   set<pair<int, int> >  col_state_var_incidence;
-                  for (set<pair<int, int> >::const_iterator row_state_var_incidence_it = row_state_var_incidence.begin(); row_state_var_incidence_it != row_state_var_incidence.end(); row_state_var_incidence_it++)
-                    col_state_var_incidence.insert(make_pair(row_state_var_incidence_it->second, row_state_var_incidence_it->first));
+                  for (const auto & row_state_var_incidence_it : row_state_var_incidence)
+                    col_state_var_incidence.insert(make_pair(row_state_var_incidence_it.second, row_state_var_incidence_it.first));
                   set<pair<int, int> >::const_iterator  col_state_var_incidence_it = col_state_var_incidence.begin();
                   diag = true;
                   int nb_diag_c = 0;
@@ -3398,8 +3383,8 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
           int size_v_index_KF = v_index_KF.size();
 
           KF_index_file.write(reinterpret_cast<char *>(&size_v_index_KF), sizeof(size_v_index_KF));
-          for (vector<index_KF>::iterator it = v_index_KF.begin(); it != v_index_KF.end(); it++)
-            KF_index_file.write(reinterpret_cast<char *>(&(*it)), sizeof(index_KF));
+          for (auto & it : v_index_KF)
+            KF_index_file.write(reinterpret_cast<char *>(&it), sizeof(index_KF));
 
           vector<index_KF> v_index_KF_2;
           int n_n_obs = n * n_obs;
@@ -3417,8 +3402,8 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
           int size_v_index_KF_2 = v_index_KF_2.size();
 
           KF_index_file.write(reinterpret_cast<char *>(&size_v_index_KF_2), sizeof(size_v_index_KF_2));
-          for (vector<index_KF>::iterator it = v_index_KF_2.begin(); it != v_index_KF_2.end(); it++)
-            KF_index_file.write(reinterpret_cast<char *>(&(*it)), sizeof(index_KF));
+          for (auto & it : v_index_KF_2)
+            KF_index_file.write(reinterpret_cast<char *>(&it), sizeof(index_KF));
           KF_index_file.close();
         }
     }
@@ -3474,24 +3459,22 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
   output << "];" << endl;
 
   // Write PacExpectationInfo
-  for (set<const PacExpectationNode *>::const_iterator it = pac_expectation_info.begin();
-       it != pac_expectation_info.end(); it++)
-    (*it)->ExprNode::writeOutput(output, oMatlabDynamicModel);
+  for (auto it : pac_expectation_info)
+    it->ExprNode::writeOutput(output, oMatlabDynamicModel);
 }
 
 map<pair<int, pair<int, int > >, expr_t>
 DynamicModel::collect_first_order_derivatives_endogenous()
 {
   map<pair<int, pair<int, int > >, expr_t> endo_derivatives;
-  for (first_derivatives_t::iterator it2 = first_derivatives.begin();
-       it2 != first_derivatives.end(); it2++)
+  for (auto & first_derivative : first_derivatives)
     {
-      if (getTypeByDerivID(it2->first.second) == eEndogenous)
+      if (getTypeByDerivID(first_derivative.first.second) == eEndogenous)
         {
-          int eq = it2->first.first;
-          int var = symbol_table.getTypeSpecificID(getSymbIDByDerivID(it2->first.second));
-          int lag = getLagByDerivID(it2->first.second);
-          endo_derivatives[make_pair(eq, make_pair(var, lag))] = it2->second;
+          int eq = first_derivative.first.first;
+          int var = symbol_table.getTypeSpecificID(getSymbIDByDerivID(first_derivative.first.second));
+          int lag = getLagByDerivID(first_derivative.first.second);
+          endo_derivatives[make_pair(eq, make_pair(var, lag))] = first_derivative.second;
         }
     }
   return endo_derivatives;
@@ -3518,12 +3501,11 @@ DynamicModel::getVarModelVariablesFromEqTags(vector<string> &var_model_eqtags,
       int eqn = -1;
       set<pair<int, int> > lhs_set, lhs_tmp_set, rhs_set;
       string eqtag (*itvareqs);
-      for (vector<pair<int, pair<string, string> > >::const_iterator iteqtag =
-             equation_tags.begin(); iteqtag != equation_tags.end(); iteqtag++)
-        if (iteqtag->second.first == "name"
-            && iteqtag->second.second == eqtag)
+      for (const auto & equation_tag : equation_tags)
+        if (equation_tag.second.first == "name"
+            && equation_tag.second.second == eqtag)
           {
-            eqn = iteqtag->first;
+            eqn = equation_tag.first;
             break;
           }
 
@@ -3534,11 +3516,10 @@ DynamicModel::getVarModelVariablesFromEqTags(vector<string> &var_model_eqtags,
         }
 
       bool nonstationary_bool = false;
-      for (vector<pair<int, pair<string, string> > >::const_iterator iteqtag =
-             equation_tags.begin(); iteqtag != equation_tags.end(); iteqtag++)
-        if (iteqtag->first == eqn)
-          if (iteqtag->second.first == "data_type"
-              && iteqtag->second.second == "nonstationary")
+      for (const auto & equation_tag : equation_tags)
+        if (equation_tag.first == eqn)
+          if (equation_tag.second.first == "data_type"
+              && equation_tag.second.second == "nonstationary")
             {
               nonstationary_bool = true;
               break;
@@ -3617,9 +3598,8 @@ DynamicModel::getVarMaxLag(StaticModel &static_model, vector<int> &eqnumber) con
     }
 
   set<expr_t> lhs_static;
-  for(set<expr_t>::const_iterator it = lhs.begin();
-      it != lhs.end(); it++)
-    lhs_static.insert((*it)->toStatic(static_model));
+  for(auto lh : lhs)
+    lhs_static.insert(lh->toStatic(static_model));
 
   int max_lag = 0;
   for (vector<int>::const_iterator it = eqnumber.begin();
@@ -3657,8 +3637,8 @@ DynamicModel::getVarLhsDiffAndInfo(vector<int> &eqnumber, vector<bool> &diff,
 void
 DynamicModel::setVarExpectationIndices(map<string, pair<SymbolList, int> > &var_model_info)
 {
-  for (size_t i = 0; i < equations.size(); i++)
-    equations[i]->setVarExpectationIndex(var_model_info);
+  for (auto & equation : equations)
+    equation->setVarExpectationIndex(var_model_info);
 }
 
 void
@@ -3668,8 +3648,8 @@ DynamicModel::addEquationsForVar(map<string, pair<SymbolList, int> > &var_model_
   map<string, int> var_endos_and_lags, model_endos_and_lags;
   for (map<string, pair<SymbolList, int> >::const_iterator it = var_model_info.begin();
        it != var_model_info.end(); it++)
-    for (size_t i = 0; i < equations.size(); i++)
-      if (equations[i]->isVarModelReferenced(it->first))
+    for (auto & equation : equations)
+      if (equation->isVarModelReferenced(it->first))
         {
           vector<string> symbol_list = it->second.first.get_symbols();
           int order = it->second.second;
@@ -3687,8 +3667,8 @@ DynamicModel::addEquationsForVar(map<string, pair<SymbolList, int> > &var_model_
     return;
 
   // Ensure that the minimum lag value exists in the model equations. If not, add an equation for it
-  for (size_t i = 0; i < equations.size(); i++)
-    equations[i]->getEndosAndMaxLags(model_endos_and_lags);
+  for (auto & equation : equations)
+    equation->getEndosAndMaxLags(model_endos_and_lags);
 
   int count = 0;
   for (map<string, int>::const_iterator it = var_endos_and_lags.begin();
@@ -3824,15 +3804,15 @@ DynamicModel::getUndiffMaxLag(StaticModel &static_model, vector<expr_t> &lhs, ve
 void
 DynamicModel::walkPacParameters()
 {
-  for (size_t i = 0; i < equations.size(); i++)
+  for (auto & equation : equations)
     {
       bool pac_encountered = false;
       pair<int, int> lhs (-1, -1);
       set<pair<int, pair<int, int> > > params_and_vars;
       set<pair<int, pair<int, int> > > ecm_params_and_vars;
-      equations[i]->walkPacParameters(pac_encountered, lhs, ecm_params_and_vars, params_and_vars);
+      equation->walkPacParameters(pac_encountered, lhs, ecm_params_and_vars, params_and_vars);
       if (pac_encountered)
-        equations[i]->addParamInfoToPac(lhs, ecm_params_and_vars, params_and_vars);
+        equation->addParamInfoToPac(lhs, ecm_params_and_vars, params_and_vars);
     }
 }
 
@@ -3851,15 +3831,14 @@ void
 DynamicModel::substitutePacExpectation()
 {
   map<const PacExpectationNode *, const BinaryOpNode *> subst_table;
-  for (map<int, expr_t>::iterator it = local_variables_table.begin();
-       it != local_variables_table.end(); it++)
-    it->second = it->second->substitutePacExpectation(subst_table);
+  for (auto & it : local_variables_table)
+    it.second = it.second->substitutePacExpectation(subst_table);
 
-  for (size_t i = 0; i < equations.size(); i++)
+  for (auto & equation : equations)
     {
-      BinaryOpNode *substeq = dynamic_cast<BinaryOpNode *>(equations[i]->substitutePacExpectation(subst_table));
+      BinaryOpNode *substeq = dynamic_cast<BinaryOpNode *>(equation->substitutePacExpectation(subst_table));
       assert(substeq != NULL);
-      equations[i] = substeq;
+      equation = substeq;
     }
 
   for (map<const PacExpectationNode *, const BinaryOpNode *>::const_iterator it = subst_table.begin();
@@ -3989,11 +3968,10 @@ void
 DynamicModel::computeXrefs()
 {
   int i = 0;
-  for (vector<BinaryOpNode *>::iterator it = equations.begin();
-       it != equations.end(); it++)
+  for (auto & equation : equations)
     {
       ExprNode::EquationInfo ei;
-      (*it)->computeXrefs(ei);
+      equation->computeXrefs(ei);
       xrefs[i++] = ei;
     }
 
@@ -4011,14 +3989,13 @@ DynamicModel::computeXrefs()
 void
 DynamicModel::computeRevXref(map<pair<int, int>, set<int> > &xrefset, const set<pair<int, int> > &eiref, int eqn)
 {
-  for (set<pair<int, int> >::const_iterator it = eiref.begin();
-       it != eiref.end(); it++)
+  for (const auto & it : eiref)
     {
       set<int> eq;
-      if (xrefset.find(*it) != xrefset.end())
-        eq = xrefset[*it];
+      if (xrefset.find(it) != xrefset.end())
+        eq = xrefset[it];
       eq.insert(eqn);
-      xrefset[*it] = eq;
+      xrefset[it] = eq;
     }
 }
 
@@ -4034,27 +4011,23 @@ DynamicModel::writeXrefs(ostream &output) const
        it != xrefs.end(); it++, i++)
     {
       output << "M_.xref1.param{" << i << "} = [ ";
-      for (set<pair<int, int> >::const_iterator it1 = it->second.param.begin();
-           it1 != it->second.param.end(); it1++)
-        output << symbol_table.getTypeSpecificID(it1->first) + 1 << " ";
+      for (const auto & it1 : it->second.param)
+        output << symbol_table.getTypeSpecificID(it1.first) + 1 << " ";
       output << "];" << endl;
 
       output << "M_.xref1.endo{" << i << "} = [ ";
-      for (set<pair<int, int> >::const_iterator it1 = it->second.endo.begin();
-           it1 != it->second.endo.end(); it1++)
-        output << "struct('id', " << symbol_table.getTypeSpecificID(it1->first) + 1 << ", 'shift', " << it1->second << ");";
+      for (const auto & it1 : it->second.endo)
+        output << "struct('id', " << symbol_table.getTypeSpecificID(it1.first) + 1 << ", 'shift', " << it1.second << ");";
       output << "];" << endl;
 
       output << "M_.xref1.exo{" << i << "} = [ ";
-      for (set<pair<int, int> >::const_iterator it1 = it->second.exo.begin();
-           it1 != it->second.exo.end(); it1++)
-        output << "struct('id', " << symbol_table.getTypeSpecificID(it1->first) + 1 << ", 'shift', " << it1->second << ");";
+      for (const auto & it1 : it->second.exo)
+        output << "struct('id', " << symbol_table.getTypeSpecificID(it1.first) + 1 << ", 'shift', " << it1.second << ");";
       output << "];" << endl;
 
       output << "M_.xref1.exo_det{" << i << "} = [ ";
-      for (set<pair<int, int> >::const_iterator it1 = it->second.exo_det.begin();
-           it1 != it->second.exo_det.end(); it1++)
-        output << "struct('id', " << symbol_table.getTypeSpecificID(it1->first) + 1 << ", 'shift', " << it1->second << ");";
+      for (const auto & it1 : it->second.exo_det)
+        output << "struct('id', " << symbol_table.getTypeSpecificID(it1.first) + 1 << ", 'shift', " << it1.second << ");";
       output << "];" << endl;
     }
 
@@ -4072,22 +4045,21 @@ void
 DynamicModel::writeRevXrefs(ostream &output, const map<pair<int, int>, set<int> > &xrefmap, const string &type) const
 {
   int last_tsid = -1;
-  for (map<pair<int, int>, set<int> >::const_iterator it = xrefmap.begin();
-       it != xrefmap.end(); it++)
+  for (const auto & it : xrefmap)
     {
-      int tsid = symbol_table.getTypeSpecificID(it->first.first) + 1;
+      int tsid = symbol_table.getTypeSpecificID(it.first.first) + 1;
       output << "M_.xref2." << type << "{" << tsid << "} = [ ";
       if (last_tsid == tsid)
         output << "M_.xref2." << type << "{" << tsid << "}; ";
       else
         last_tsid = tsid;
 
-      for (set<int>::const_iterator it1 = it->second.begin();
-           it1 != it->second.end(); it1++)
+      for (set<int>::const_iterator it1 = it.second.begin();
+           it1 != it.second.end(); it1++)
         if (type == "param")
           output << *it1 + 1 << " ";
         else
-          output << "struct('shift', " << it->first.second << ", 'eq', " << *it1+1 << ");";
+          output << "struct('shift', " << it.first.second << ", 'eq', " << *it1+1 << ");";
       output << "];" << endl;
     }
 }
@@ -4236,17 +4208,16 @@ DynamicModel::collect_block_first_order_derivatives()
   exo_max_leadlag_block = vector<pair<int, int> >(nb_blocks, make_pair(0, 0));
   exo_det_max_leadlag_block = vector<pair<int, int> >(nb_blocks, make_pair(0, 0));
   max_leadlag_block = vector<pair<int, int> >(nb_blocks, make_pair(0, 0));
-  for (first_derivatives_t::iterator it2 = first_derivatives.begin();
-       it2 != first_derivatives.end(); it2++)
+  for (auto & first_derivative : first_derivatives)
     {
-      int eq = it2->first.first;
-      int var = symbol_table.getTypeSpecificID(getSymbIDByDerivID(it2->first.second));
-      int lag = getLagByDerivID(it2->first.second);
+      int eq = first_derivative.first.first;
+      int var = symbol_table.getTypeSpecificID(getSymbIDByDerivID(first_derivative.first.second));
+      int lag = getLagByDerivID(first_derivative.first.second);
       int block_eq = equation_2_block[eq];
       int block_var = 0;
       derivative_t tmp_derivative;
       lag_var_t lag_var;
-      switch (getTypeByDerivID(it2->first.second))
+      switch (getTypeByDerivID(first_derivative.first.second))
         {
         case eEndogenous:
           block_var = variable_2_block[var];
@@ -4369,9 +4340,8 @@ DynamicModel::collectBlockVariables()
       for (lag_var_t::const_iterator it = exo_block[block].begin(); it != exo_block[block].end(); it++)
         {
           int lag = it->first;
-          for (var_t::const_iterator it2 = it->second.begin(); it2 != it->second.end(); it2++)
+          for (int var : it->second)
             {
-              int var = *it2;
               tmp_var_exo.insert(var);
               if (prev_var != var || prev_lag != lag)
                 {
@@ -4457,15 +4427,15 @@ DynamicModel::writeAuxVarRecursiveDefinitions(ostream &output, ExprNodeOutputTyp
   deriv_node_temp_terms_t tef_terms;
   temporary_terms_t temporary_terms;
   temporary_terms_idxs_t temporary_terms_idxs;
-  for (int i = 0; i < (int) aux_equations.size(); i++)
-    if (dynamic_cast<ExprNode *>(aux_equations[i])->containsExternalFunction())
-      dynamic_cast<ExprNode *>(aux_equations[i])->writeExternalFunctionOutput(output, output_type,
+  for (auto aux_equation : aux_equations)
+    if (dynamic_cast<ExprNode *>(aux_equation)->containsExternalFunction())
+      dynamic_cast<ExprNode *>(aux_equation)->writeExternalFunctionOutput(output, output_type,
                                                                               temporary_terms,
                                                                               temporary_terms_idxs,
                                                                               tef_terms);
-  for (int i = 0; i < (int) aux_equations.size(); i++)
+  for (auto aux_equation : aux_equations)
     {
-      dynamic_cast<ExprNode *>(aux_equations[i])->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
+      dynamic_cast<ExprNode *>(aux_equation)->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
       output << ";" << endl;
     }
 }
@@ -4493,25 +4463,22 @@ DynamicModel::cloneDynamic(DynamicModel &dynamic_model) const
   assert(&symbol_table == &dynamic_model.symbol_table);
 
   // Convert model local variables (need to be done first)
-  for (vector<int>::const_iterator it = local_variables_vector.begin();
-       it != local_variables_vector.end(); it++)
-    dynamic_model.AddLocalVariable(*it, local_variables_table.find(*it)->second->cloneDynamic(dynamic_model));
+  for (int it : local_variables_vector)
+    dynamic_model.AddLocalVariable(it, local_variables_table.find(it)->second->cloneDynamic(dynamic_model));
 
   // Convert equations
   for (size_t i = 0; i < equations.size(); i++)
     {
       vector<pair<string, string> > eq_tags;
-      for (vector<pair<int, pair<string, string> > >::const_iterator it = equation_tags.begin();
-           it != equation_tags.end(); ++it)
-        if (it->first == (int)i)
-          eq_tags.push_back(it->second);
+      for (const auto & equation_tag : equation_tags)
+        if (equation_tag.first == (int)i)
+          eq_tags.push_back(equation_tag.second);
       dynamic_model.addEquation(equations[i]->cloneDynamic(dynamic_model), equations_lineno[i], eq_tags);
     }
 
   // Convert auxiliary equations
-  for (deque<BinaryOpNode *>::const_iterator it = aux_equations.begin();
-       it != aux_equations.end(); it++)
-    dynamic_model.addAuxEquation((*it)->cloneDynamic(dynamic_model));
+  for (auto aux_equation : aux_equations)
+    dynamic_model.addAuxEquation(aux_equation->cloneDynamic(dynamic_model));
 
   // Convert static_only equations
   for (size_t i = 0; i < static_only_equations.size(); i++)
@@ -4554,13 +4521,12 @@ DynamicModel::computeRamseyPolicyFOCs(const StaticModel &static_model, const boo
   set<pair<int, int> > dynvars;
   int max_eq_lead = 0;
   int max_eq_lag = 0;
-  for (int i = 0; i < (int) equations.size(); i++)
-    equations[i]->collectDynamicVariables(eEndogenous, dynvars);
+  for (auto & equation : equations)
+    equation->collectDynamicVariables(eEndogenous, dynvars);
 
-  for (set<pair<int, int> >::const_iterator it = dynvars.begin();
-       it != dynvars.end(); it++)
+  for (const auto & dynvar : dynvars)
     {
-      int lag = it->second;
+      int lag = dynvar.second;
       if (max_eq_lead < lag)
         max_eq_lead = lag;
       else if (-max_eq_lag > lag)
@@ -4606,8 +4572,8 @@ DynamicModel::computeRamseyPolicyFOCs(const StaticModel &static_model, const boo
 
   // Add new equations
   equations.clear();
-  for (int i = 0; i < (int) neweqs.size(); i++)
-    addEquation(neweqs[i], -1);
+  for (auto & neweq : neweqs)
+    addEquation(neweq, -1);
 }
 
 void
@@ -4618,9 +4584,8 @@ DynamicModel::toStatic(StaticModel &static_model) const
   assert(&symbol_table == &static_model.symbol_table);
 
   // Convert model local variables (need to be done first)
-  for (vector<int>::const_iterator it = local_variables_vector.begin();
-       it != local_variables_vector.end(); it++)
-    static_model.AddLocalVariable(*it, local_variables_table.find(*it)->second->toStatic(static_model));
+  for (int it : local_variables_vector)
+    static_model.AddLocalVariable(it, local_variables_table.find(it)->second->toStatic(static_model));
 
   // Convert equations
   int static_only_index = 0;
@@ -4629,12 +4594,11 @@ DynamicModel::toStatic(StaticModel &static_model) const
       // Detect if equation is marked [dynamic]
       bool is_dynamic_only = false;
       vector<pair<string, string> > eq_tags;
-      for (vector<pair<int, pair<string, string> > >::const_iterator it = equation_tags.begin();
-           it != equation_tags.end(); ++it)
-        if (it->first == i)
+      for (const auto & equation_tag : equation_tags)
+        if (equation_tag.first == i)
           {
-            eq_tags.push_back(it->second);
-            if (it->second.first == "dynamic")
+            eq_tags.push_back(equation_tag.second);
+            if (equation_tag.second.first == "dynamic")
               is_dynamic_only = true;
           }
 
@@ -4657,9 +4621,8 @@ DynamicModel::toStatic(StaticModel &static_model) const
     }
 
   // Convert auxiliary equations
-  for (deque<BinaryOpNode *>::const_iterator it = aux_equations.begin();
-       it != aux_equations.end(); it++)
-    static_model.addAuxEquation((*it)->toStatic(static_model));
+  for (auto aux_equation : aux_equations)
+    static_model.addAuxEquation(aux_equation->toStatic(static_model));
 }
 
 bool
@@ -4672,8 +4635,8 @@ set<int>
 DynamicModel::findUnusedEndogenous()
 {
   set<int> usedEndo, unusedEndo;
-  for (int i = 0; i < (int) equations.size(); i++)
-    equations[i]->collectVariables(eEndogenous, usedEndo);
+  for (auto & equation : equations)
+    equation->collectVariables(eEndogenous, usedEndo);
   set<int> allEndo = symbol_table.getEndogenous();
   set_difference(allEndo.begin(), allEndo.end(),
                  usedEndo.begin(), usedEndo.end(),
@@ -4685,8 +4648,8 @@ set<int>
 DynamicModel::findUnusedExogenous()
 {
   set<int> usedExo, unusedExo, unobservedExo;
-  for (int i = 0; i < (int) equations.size(); i++)
-    equations[i]->collectVariables(eExogenous, usedExo);
+  for (auto & equation : equations)
+    equation->collectVariables(eExogenous, usedExo);
   set<int> observedExo = symbol_table.getObservedExogenous();
   set<int> allExo = symbol_table.getExogenous();
   set_difference(allExo.begin(), allExo.end(),
@@ -4703,18 +4666,17 @@ DynamicModel::setLeadsLagsOrig()
 {
   set<pair<int, int> > dynvars;
 
-  for (int i = 0; i < (int) equations.size(); i++)
+  for (auto & equation : equations)
     {
-      equations[i]->collectDynamicVariables(eEndogenous, dynvars);
-      equations[i]->collectDynamicVariables(eExogenous, dynvars);
-      equations[i]->collectDynamicVariables(eExogenousDet, dynvars);
+      equation->collectDynamicVariables(eEndogenous, dynvars);
+      equation->collectDynamicVariables(eExogenous, dynvars);
+      equation->collectDynamicVariables(eExogenousDet, dynvars);
     }
 
-    for (set<pair<int, int> >::const_iterator it = dynvars.begin();
-         it != dynvars.end(); it++)
+    for (const auto & dynvar : dynvars)
     {
-      int lag = it->second;
-      SymbolType type = symbol_table.getType(it->first);
+      int lag = dynvar.second;
+      SymbolType type = symbol_table.getType(dynvar.first);
 
       if (max_lead_orig < lag)
         max_lead_orig= lag;
@@ -4752,25 +4714,24 @@ DynamicModel::computeDerivIDs()
 {
   set<pair<int, int> > dynvars;
 
-  for (int i = 0; i < (int) equations.size(); i++)
-    equations[i]->collectDynamicVariables(eEndogenous, dynvars);
+  for (auto & equation : equations)
+    equation->collectDynamicVariables(eEndogenous, dynvars);
 
   dynJacobianColsNbr = dynvars.size();
 
-  for (int i = 0; i < (int) equations.size(); i++)
+  for (auto & equation : equations)
     {
-      equations[i]->collectDynamicVariables(eExogenous, dynvars);
-      equations[i]->collectDynamicVariables(eExogenousDet, dynvars);
-      equations[i]->collectDynamicVariables(eParameter, dynvars);
-      equations[i]->collectDynamicVariables(eTrend, dynvars);
-      equations[i]->collectDynamicVariables(eLogTrend, dynvars);
+      equation->collectDynamicVariables(eExogenous, dynvars);
+      equation->collectDynamicVariables(eExogenousDet, dynvars);
+      equation->collectDynamicVariables(eParameter, dynvars);
+      equation->collectDynamicVariables(eTrend, dynvars);
+      equation->collectDynamicVariables(eLogTrend, dynvars);
     }
 
-  for (set<pair<int, int> >::const_iterator it = dynvars.begin();
-       it != dynvars.end(); it++)
+  for (const auto & dynvar : dynvars)
     {
-      int lag = it->second;
-      SymbolType type = symbol_table.getType(it->first);
+      int lag = dynvar.second;
+      SymbolType type = symbol_table.getType(dynvar.first);
 
       /* Setting maximum and minimum lags.
 
@@ -4809,8 +4770,8 @@ DynamicModel::computeDerivIDs()
       // Create a new deriv_id
       int deriv_id = deriv_id_table.size();
 
-      deriv_id_table[*it] = deriv_id;
-      inv_deriv_id_table.push_back(*it);
+      deriv_id_table[dynvar] = deriv_id;
+      inv_deriv_id_table.push_back(dynvar);
     }
 }
 
@@ -4978,12 +4939,11 @@ DynamicModel::writeParamsDerivativesFile(const string &basename, bool julia) con
 
   writeTemporaryTerms(params_derivs_temporary_terms, {}, params_derivs_temporary_terms_idxs, model_output, output_type, tef_terms);
 
-  for (first_derivatives_t::const_iterator it = residuals_params_derivatives.begin();
-       it != residuals_params_derivatives.end(); it++)
+  for (const auto & residuals_params_derivative : residuals_params_derivatives)
     {
-      int eq = it->first.first;
-      int param = it->first.second;
-      expr_t d1 = it->second;
+      int eq = residuals_params_derivative.first.first;
+      int param = residuals_params_derivative.first.second;
+      expr_t d1 = residuals_params_derivative.second;
 
       int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1;
 
@@ -4993,13 +4953,12 @@ DynamicModel::writeParamsDerivativesFile(const string &basename, bool julia) con
       jacobian_output << ";" << endl;
     }
 
-  for (second_derivatives_t::const_iterator it = jacobian_params_derivatives.begin();
-       it != jacobian_params_derivatives.end(); it++)
+  for (const auto & jacobian_params_derivative : jacobian_params_derivatives)
     {
-      int eq = it->first.first;
-      int var = it->first.second.first;
-      int param = it->first.second.second;
-      expr_t d2 = it->second;
+      int eq = jacobian_params_derivative.first.first;
+      int var = jacobian_params_derivative.first.second.first;
+      int param = jacobian_params_derivative.first.second.second;
+      expr_t d2 = jacobian_params_derivative.second;
 
       int var_col = getDynJacobianCol(var) + 1;
       int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1;
@@ -5248,13 +5207,12 @@ DynamicModel::substituteLeadLagInternal(aux_var_t type, bool deterministic_model
 
   // Substitute in used model local variables
   set<int> used_local_vars;
-  for (size_t i = 0; i < equations.size(); i++)
-    equations[i]->collectVariables(eModelLocalVariable, used_local_vars);
+  for (auto & equation : equations)
+    equation->collectVariables(eModelLocalVariable, used_local_vars);
 
-  for (set<int>::const_iterator it = used_local_vars.begin();
-       it != used_local_vars.end(); ++it)
+  for (int used_local_var : used_local_vars)
     {
-      const expr_t value = local_variables_table.find(*it)->second;
+      const expr_t value = local_variables_table.find(used_local_var)->second;
       expr_t subst;
       switch (type)
         {
@@ -5277,29 +5235,29 @@ DynamicModel::substituteLeadLagInternal(aux_var_t type, bool deterministic_model
           cerr << "DynamicModel::substituteLeadLagInternal: impossible case" << endl;
           exit(EXIT_FAILURE);
         }
-      local_variables_table[*it] = subst;
+      local_variables_table[used_local_var] = subst;
     }
 
   // Substitute in equations
-  for (int i = 0; i < (int) equations.size(); i++)
+  for (auto & equation : equations)
     {
       expr_t subst;
       switch (type)
         {
         case avEndoLead:
-          subst = equations[i]->substituteEndoLeadGreaterThanTwo(subst_table, neweqs, deterministic_model);
+          subst = equation->substituteEndoLeadGreaterThanTwo(subst_table, neweqs, deterministic_model);
           break;
         case avEndoLag:
-          subst = equations[i]->substituteEndoLagGreaterThanTwo(subst_table, neweqs);
+          subst = equation->substituteEndoLagGreaterThanTwo(subst_table, neweqs);
           break;
         case avExoLead:
-          subst = equations[i]->substituteExoLead(subst_table, neweqs, deterministic_model);
+          subst = equation->substituteExoLead(subst_table, neweqs, deterministic_model);
           break;
         case avExoLag:
-          subst = equations[i]->substituteExoLag(subst_table, neweqs);
+          subst = equation->substituteExoLag(subst_table, neweqs);
           break;
         case avDiffForward:
-          subst = equations[i]->differentiateForwardVars(subset, subst_table, neweqs);
+          subst = equation->differentiateForwardVars(subset, subst_table, neweqs);
           break;
         default:
           cerr << "DynamicModel::substituteLeadLagInternal: impossible case" << endl;
@@ -5307,32 +5265,32 @@ DynamicModel::substituteLeadLagInternal(aux_var_t type, bool deterministic_model
         }
       BinaryOpNode *substeq = dynamic_cast<BinaryOpNode *>(subst);
       assert(substeq != NULL);
-      equations[i] = substeq;
+      equation = substeq;
     }
 
   // Substitute in aux_equations
   // Without this loop, the auxiliary equations in equations
   // will diverge from those in aux_equations
-  for (int i = 0; i < (int) aux_equations.size(); i++)
+  for (auto & aux_equation : aux_equations)
     {
       expr_t subst;
       switch (type)
         {
         case avEndoLead:
-          subst = aux_equations[i]->substituteEndoLeadGreaterThanTwo(subst_table,
+          subst = aux_equation->substituteEndoLeadGreaterThanTwo(subst_table,
                                                                      neweqs, deterministic_model);
           break;
         case avEndoLag:
-          subst = aux_equations[i]->substituteEndoLagGreaterThanTwo(subst_table, neweqs);
+          subst = aux_equation->substituteEndoLagGreaterThanTwo(subst_table, neweqs);
           break;
         case avExoLead:
-          subst = aux_equations[i]->substituteExoLead(subst_table, neweqs, deterministic_model);
+          subst = aux_equation->substituteExoLead(subst_table, neweqs, deterministic_model);
           break;
         case avExoLag:
-          subst = aux_equations[i]->substituteExoLag(subst_table, neweqs);
+          subst = aux_equation->substituteExoLag(subst_table, neweqs);
           break;
         case avDiffForward:
-          subst = aux_equations[i]->differentiateForwardVars(subset, subst_table, neweqs);
+          subst = aux_equation->differentiateForwardVars(subset, subst_table, neweqs);
           break;
         default:
           cerr << "DynamicModel::substituteLeadLagInternal: impossible case" << endl;
@@ -5340,32 +5298,32 @@ DynamicModel::substituteLeadLagInternal(aux_var_t type, bool deterministic_model
         }
       BinaryOpNode *substeq = dynamic_cast<BinaryOpNode *>(subst);
       assert(substeq != NULL);
-      aux_equations[i] = substeq;
+      aux_equation = substeq;
     }
 
  // Substitute in diff_aux_equations
   // Without this loop, the auxiliary equations in equations
   // will diverge from those in diff_aux_equations
-  for (int i = 0; i < (int) diff_aux_equations.size(); i++)
+  for (auto & diff_aux_equation : diff_aux_equations)
     {
       expr_t subst;
       switch (type)
         {
         case avEndoLead:
-          subst = diff_aux_equations[i]->substituteEndoLeadGreaterThanTwo(subst_table,
+          subst = diff_aux_equation->substituteEndoLeadGreaterThanTwo(subst_table,
                                                                      neweqs, deterministic_model);
           break;
         case avEndoLag:
-          subst = diff_aux_equations[i]->substituteEndoLagGreaterThanTwo(subst_table, neweqs);
+          subst = diff_aux_equation->substituteEndoLagGreaterThanTwo(subst_table, neweqs);
           break;
         case avExoLead:
-          subst = diff_aux_equations[i]->substituteExoLead(subst_table, neweqs, deterministic_model);
+          subst = diff_aux_equation->substituteExoLead(subst_table, neweqs, deterministic_model);
           break;
         case avExoLag:
-          subst = diff_aux_equations[i]->substituteExoLag(subst_table, neweqs);
+          subst = diff_aux_equation->substituteExoLag(subst_table, neweqs);
           break;
         case avDiffForward:
-          subst = diff_aux_equations[i]->differentiateForwardVars(subset, subst_table, neweqs);
+          subst = diff_aux_equation->differentiateForwardVars(subset, subst_table, neweqs);
           break;
         default:
           cerr << "DynamicModel::substituteLeadLagInternal: impossible case" << endl;
@@ -5373,12 +5331,12 @@ DynamicModel::substituteLeadLagInternal(aux_var_t type, bool deterministic_model
         }
       BinaryOpNode *substeq = dynamic_cast<BinaryOpNode *>(subst);
       assert(substeq != NULL);
-      diff_aux_equations[i] = substeq;
+      diff_aux_equation = substeq;
     }
 
   // Add new equations
-  for (int i = 0; i < (int) neweqs.size(); i++)
-    addEquation(neweqs[i], -1);
+  for (auto & neweq : neweqs)
+    addEquation(neweq, -1);
 
   // Order of auxiliary variable definition equations:
   //  - expectation (entered before this function is called)
@@ -5420,8 +5378,8 @@ DynamicModel::substituteLeadLagInternal(aux_var_t type, bool deterministic_model
 void
 DynamicModel::substituteAdl()
 {
-  for (int i = 0; i < (int) equations.size(); i++)
-    equations[i] = dynamic_cast<BinaryOpNode *>(equations[i]->substituteAdl());
+  for (auto & equation : equations)
+    equation = dynamic_cast<BinaryOpNode *>(equation->substituteAdl());
 }
 
 void
@@ -5429,32 +5387,30 @@ DynamicModel::substituteUnaryOps(StaticModel &static_model)
 {
   diff_table_t nodes;
   // Find matching unary ops that may be outside of diffs (i.e., those with different lags)
-  for (map<int, expr_t>::iterator it = local_variables_table.begin();
-       it != local_variables_table.end(); it++)
-    it->second->findUnaryOpNodesForAuxVarCreation(static_model, nodes);
+  for (auto & it : local_variables_table)
+    it.second->findUnaryOpNodesForAuxVarCreation(static_model, nodes);
 
-  for (int i = 0; i < (int) equations.size(); i++)
-    equations[i]->findUnaryOpNodesForAuxVarCreation(static_model, nodes);
+  for (auto & equation : equations)
+    equation->findUnaryOpNodesForAuxVarCreation(static_model, nodes);
 
   // Substitute in model local variables
   ExprNode::subst_table_t subst_table;
   vector<BinaryOpNode *> neweqs;
-  for (map<int, expr_t>::iterator it = local_variables_table.begin();
-       it != local_variables_table.end(); it++)
-    it->second = it->second->substituteUnaryOpNodes(static_model, nodes, subst_table, neweqs);
+  for (auto & it : local_variables_table)
+    it.second = it.second->substituteUnaryOpNodes(static_model, nodes, subst_table, neweqs);
 
   // Substitute in equations
-  for (int i = 0; i < (int) equations.size(); i++)
+  for (auto & equation : equations)
     {
-      BinaryOpNode *substeq = dynamic_cast<BinaryOpNode *>(equations[i]->
+      BinaryOpNode *substeq = dynamic_cast<BinaryOpNode *>(equation->
                                                            substituteUnaryOpNodes(static_model, nodes, subst_table, neweqs));
       assert(substeq != NULL);
-      equations[i] = substeq;
+      equation = substeq;
     }
 
   // Add new equations
-  for (int i = 0; i < (int) neweqs.size(); i++)
-    addEquation(neweqs[i], -1);
+  for (auto & neweq : neweqs)
+    addEquation(neweq, -1);
 
   copy(neweqs.begin(), neweqs.end(), back_inserter(diff_aux_equations));
 
@@ -5467,31 +5423,29 @@ DynamicModel::substituteDiff(StaticModel &static_model, ExprNode::subst_table_t
 {
   // Find diff Nodes
   diff_table_t diff_table;
-  for (map<int, expr_t>::iterator it = local_variables_table.begin();
-       it != local_variables_table.end(); it++)
-    it->second->findDiffNodes(static_model, diff_table);
+  for (auto & it : local_variables_table)
+    it.second->findDiffNodes(static_model, diff_table);
 
-  for (int i = 0; i < (int) equations.size(); i++)
-    equations[i]->findDiffNodes(static_model, diff_table);
+  for (auto & equation : equations)
+    equation->findDiffNodes(static_model, diff_table);
 
   // Substitute in model local variables
   vector<BinaryOpNode *> neweqs;
-  for (map<int, expr_t>::iterator it = local_variables_table.begin();
-       it != local_variables_table.end(); it++)
-    it->second = it->second->substituteDiff(static_model, diff_table, diff_subst_table, neweqs);
+  for (auto & it : local_variables_table)
+    it.second = it.second->substituteDiff(static_model, diff_table, diff_subst_table, neweqs);
 
   // Substitute in equations
-  for (int i = 0; i < (int) equations.size(); i++)
+  for (auto & equation : equations)
     {
-      BinaryOpNode *substeq = dynamic_cast<BinaryOpNode *>(equations[i]->
+      BinaryOpNode *substeq = dynamic_cast<BinaryOpNode *>(equation->
                                                            substituteDiff(static_model, diff_table, diff_subst_table, neweqs));
       assert(substeq != NULL);
-      equations[i] = substeq;
+      equation = substeq;
     }
 
   // Add new equations
-  for (int i = 0; i < (int) neweqs.size(); i++)
-    addEquation(neweqs[i], -1);
+  for (auto & neweq : neweqs)
+    addEquation(neweq, -1);
 
   copy(neweqs.begin(), neweqs.end(), back_inserter(diff_aux_equations));
 
@@ -5512,21 +5466,20 @@ DynamicModel::substituteExpectation(bool partial_information_model)
   vector<BinaryOpNode *> neweqs;
 
   // Substitute in model local variables
-  for (map<int, expr_t>::iterator it = local_variables_table.begin();
-       it != local_variables_table.end(); it++)
-    it->second = it->second->substituteExpectation(subst_table, neweqs, partial_information_model);
+  for (auto & it : local_variables_table)
+    it.second = it.second->substituteExpectation(subst_table, neweqs, partial_information_model);
 
   // Substitute in equations
-  for (int i = 0; i < (int) equations.size(); i++)
+  for (auto & equation : equations)
     {
-      BinaryOpNode *substeq = dynamic_cast<BinaryOpNode *>(equations[i]->substituteExpectation(subst_table, neweqs, partial_information_model));
+      BinaryOpNode *substeq = dynamic_cast<BinaryOpNode *>(equation->substituteExpectation(subst_table, neweqs, partial_information_model));
       assert(substeq != NULL);
-      equations[i] = substeq;
+      equation = substeq;
     }
 
   // Add new equations
-  for (int i = 0; i < (int) neweqs.size(); i++)
-    addEquation(neweqs[i], -1);
+  for (auto & neweq : neweqs)
+    addEquation(neweq, -1);
 
   // Add the new set of equations at the *beginning* of aux_equations
   copy(neweqs.rbegin(), neweqs.rend(), front_inserter(aux_equations));
@@ -5543,15 +5496,14 @@ DynamicModel::substituteExpectation(bool partial_information_model)
 void
 DynamicModel::transformPredeterminedVariables()
 {
-  for (map<int, expr_t>::iterator it = local_variables_table.begin();
-       it != local_variables_table.end(); it++)
-    it->second = it->second->decreaseLeadsLagsPredeterminedVariables();
+  for (auto & it : local_variables_table)
+    it.second = it.second->decreaseLeadsLagsPredeterminedVariables();
 
-  for (int i = 0; i < (int) equations.size(); i++)
+  for (auto & equation : equations)
     {
-      BinaryOpNode *substeq = dynamic_cast<BinaryOpNode *>(equations[i]->decreaseLeadsLagsPredeterminedVariables());
+      BinaryOpNode *substeq = dynamic_cast<BinaryOpNode *>(equation->decreaseLeadsLagsPredeterminedVariables());
       assert(substeq != NULL);
-      equations[i] = substeq;
+      equation = substeq;
     }
 }
 
@@ -5561,29 +5513,29 @@ DynamicModel::detrendEquations()
   // We go backwards in the list of trend_vars, to deal correctly with I(2) processes
   for (nonstationary_symbols_map_t::const_reverse_iterator it = nonstationary_symbols_map.rbegin();
        it != nonstationary_symbols_map.rend(); ++it)
-    for (int i = 0; i < (int) equations.size(); i++)
+    for (auto & equation : equations)
       {
-        BinaryOpNode *substeq = dynamic_cast<BinaryOpNode *>(equations[i]->detrend(it->first, it->second.first, it->second.second));
+        BinaryOpNode *substeq = dynamic_cast<BinaryOpNode *>(equation->detrend(it->first, it->second.first, it->second.second));
         assert(substeq != NULL);
-        equations[i] = dynamic_cast<BinaryOpNode *>(substeq);
+        equation = dynamic_cast<BinaryOpNode *>(substeq);
       }
 
-  for (int i = 0; i < (int) equations.size(); i++)
+  for (auto & equation : equations)
     {
-      BinaryOpNode *substeq = dynamic_cast<BinaryOpNode *>(equations[i]->removeTrendLeadLag(trend_symbols_map));
+      BinaryOpNode *substeq = dynamic_cast<BinaryOpNode *>(equation->removeTrendLeadLag(trend_symbols_map));
       assert(substeq != NULL);
-      equations[i] = dynamic_cast<BinaryOpNode *>(substeq);
+      equation = dynamic_cast<BinaryOpNode *>(substeq);
     }
 }
 
 void
 DynamicModel::removeTrendVariableFromEquations()
 {
-  for (int i = 0; i < (int) equations.size(); i++)
+  for (auto & equation : equations)
     {
-      BinaryOpNode *substeq = dynamic_cast<BinaryOpNode *>(equations[i]->replaceTrendVar());
+      BinaryOpNode *substeq = dynamic_cast<BinaryOpNode *>(equation->replaceTrendVar());
       assert(substeq != NULL);
-      equations[i] = dynamic_cast<BinaryOpNode *>(substeq);
+      equation = dynamic_cast<BinaryOpNode *>(substeq);
     }
 }
 
@@ -5597,15 +5549,14 @@ void
 DynamicModel::fillEvalContext(eval_context_t &eval_context) const
 {
   // First, auxiliary variables
-  for (deque<BinaryOpNode *>::const_iterator it = aux_equations.begin();
-       it != aux_equations.end(); it++)
+  for (auto aux_equation : aux_equations)
     {
-      assert((*it)->get_op_code() == oEqual);
-      VariableNode *auxvar = dynamic_cast<VariableNode *>((*it)->get_arg1());
+      assert(aux_equation->get_op_code() == oEqual);
+      VariableNode *auxvar = dynamic_cast<VariableNode *>(aux_equation->get_arg1());
       assert(auxvar != NULL);
       try
         {
-          double val = (*it)->get_arg2()->eval(eval_context);
+          double val = aux_equation->get_arg2()->eval(eval_context);
           eval_context[auxvar->get_symb_id()] = val;
         }
       catch (ExprNode::EvalException &e)
@@ -5615,14 +5566,13 @@ DynamicModel::fillEvalContext(eval_context_t &eval_context) const
     }
 
   // Second, model local variables
-  for (map<int, expr_t>::const_iterator it = local_variables_table.begin();
-       it != local_variables_table.end(); it++)
+  for (auto it : local_variables_table)
     {
       try
         {
-          const expr_t expression = it->second;
+          const expr_t expression = it.second;
           double val = expression->eval(eval_context);
-          eval_context[it->first] = val;
+          eval_context[it.first] = val;
         }
       catch (ExprNode::EvalException &e)
         {
@@ -5657,8 +5607,8 @@ DynamicModel::addStaticOnlyEquation(expr_t eq, int lineno, const vector<pair<str
   assert(beq != NULL && beq->get_op_code() == oEqual);
 
   vector<pair<string, string> > soe_eq_tags;
-  for (size_t i = 0; i < eq_tags.size(); i++)
-    soe_eq_tags.push_back(eq_tags[i]);
+  for (const auto & eq_tag : eq_tags)
+    soe_eq_tags.push_back(eq_tag);
 
   static_only_equations.push_back(beq);
   static_only_equations_lineno.push_back(lineno);
@@ -5676,10 +5626,9 @@ DynamicModel::dynamicOnlyEquationsNbr() const
 {
   set<int> eqs;
 
-  for (vector<pair<int, pair<string, string> > >::const_iterator it = equation_tags.begin();
-       it != equation_tags.end(); ++it)
-    if (it->second.first == "dynamic")
-      eqs.insert(it->first);
+  for (const auto & equation_tag : equation_tags)
+    if (equation_tag.second.first == "dynamic")
+      eqs.insert(equation_tag.first);
 
   return eqs.size();
 }
@@ -5696,10 +5645,10 @@ DynamicModel::isChecksumMatching(const string &basename) const
   std::stringstream buffer;
 
   // Write equation tags
-  for (size_t i = 0; i < equation_tags.size(); i++)
-    buffer << "  " << equation_tags[i].first + 1
-           << equation_tags[i].second.first
-           << equation_tags[i].second.second;
+  for (const auto & equation_tag : equation_tags)
+    buffer << "  " << equation_tag.first + 1
+           << equation_tag.second.first
+           << equation_tag.second.second;
 
   ExprNodeOutputType buffer_type = oCDynamicModel;
 
@@ -5982,12 +5931,11 @@ DynamicModel::writeFirstDerivativesC(const string &basename, bool cuda) const
                     << "{" << endl;
 
   // Writing Jacobian
-  for (first_derivatives_t::const_iterator it = first_derivatives.begin();
-       it != first_derivatives.end(); it++)
+  for (const auto & first_derivative : first_derivatives)
     {
-      int eq = it->first.first;
-      int var = it->first.second;
-      expr_t d1 = it->second;
+      int eq = first_derivative.first.first;
+      int var = first_derivative.first.second;
+      expr_t d1 = first_derivative.second;
 
       jacobianHelper(mDynamicModelFile, eq, getDynJacobianCol(var), oCDynamicModel);
       mDynamicModelFile << "=";
@@ -6045,11 +5993,10 @@ DynamicModel::writeFirstDerivativesC_csr(const string &basename, bool cuda) cons
 
   // Indexing derivatives in column order
   vector<derivative> D;
-  for (first_derivatives_t::const_iterator it = first_derivatives.begin();
-       it != first_derivatives.end(); it++)
+  for (const auto & first_derivative : first_derivatives)
     {
-      int eq = it->first.first;
-      int dynvar = it->first.second;
+      int eq = first_derivative.first.first;
+      int dynvar = first_derivative.first.second;
       int lag = getLagByDerivID(dynvar);
       int symb_id = getSymbIDByDerivID(dynvar);
       SymbolType type = getTypeByDerivID(dynvar);
@@ -6070,7 +6017,7 @@ DynamicModel::writeFirstDerivativesC_csr(const string &basename, bool cuda) cons
           std::cerr << "This case shouldn't happen" << std::endl;
           exit(EXIT_FAILURE);
         }
-      derivative deriv(col_id + eq *cols_nbr, col_id, eq, it->second);
+      derivative deriv(col_id + eq *cols_nbr, col_id, eq, first_derivative.second);
       D.push_back(deriv);
     }
   sort(D.begin(), D.end(), derivative_less_than());
@@ -6094,9 +6041,9 @@ DynamicModel::writeFirstDerivativesC_csr(const string &basename, bool cuda) cons
   // row_ptr must point to the relative address of the first element of the row
   int cumsum = 0;
   mDynamicModelFile << "int row_ptr_data[" <<  row_ptr.size() + 1 << "] = { 0";
-  for (vector<int>::iterator it = row_ptr.begin(); it != row_ptr.end(); ++it)
+  for (int & it : row_ptr)
     {
-      cumsum += *it;
+      cumsum += it;
       mDynamicModelFile << ", " << cumsum;
     }
   mDynamicModelFile << "};" << endl
@@ -6148,24 +6095,23 @@ DynamicModel::writeSecondDerivativesC_csr(const string &basename, bool cuda) con
   // Indexing derivatives in column order
   vector<derivative> D;
   int hessianColsNbr = dynJacobianColsNbr*dynJacobianColsNbr;
-  for (second_derivatives_t::const_iterator it = second_derivatives.begin();
-       it != second_derivatives.end(); it++)
+  for (const auto & second_derivative : second_derivatives)
     {
-      int eq = it->first.first;
-      int var1 = it->first.second.first;
-      int var2 = it->first.second.second;
+      int eq = second_derivative.first.first;
+      int var1 = second_derivative.first.second.first;
+      int var2 = second_derivative.first.second.second;
 
       int id1 = getDynJacobianCol(var1);
       int id2 = getDynJacobianCol(var2);
 
       int col_nb = id1 * dynJacobianColsNbr + id2;
 
-      derivative deriv(col_nb + eq *hessianColsNbr, col_nb, eq, it->second);
+      derivative deriv(col_nb + eq *hessianColsNbr, col_nb, eq, second_derivative.second);
       D.push_back(deriv);
       if (id1 != id2)
         {
           col_nb = id2 * dynJacobianColsNbr + id1;
-          derivative deriv(col_nb + eq *hessianColsNbr, col_nb, eq, it->second);
+          derivative deriv(col_nb + eq *hessianColsNbr, col_nb, eq, second_derivative.second);
           D.push_back(deriv);
         }
     }
@@ -6190,9 +6136,9 @@ DynamicModel::writeSecondDerivativesC_csr(const string &basename, bool cuda) con
   // row_ptr must point to the relative address of the first element of the row
   int cumsum = 0;
   mDynamicModelFile << "row_ptr = [ 0";
-  for (vector<int>::iterator it = row_ptr.begin(); it != row_ptr.end(); ++it)
+  for (int & it : row_ptr)
     {
-      cumsum += *it;
+      cumsum += it;
       mDynamicModelFile << ", " << cumsum;
     }
   mDynamicModelFile << "];" << endl;
@@ -6244,13 +6190,12 @@ DynamicModel::writeThirdDerivativesC_csr(const string &basename, bool cuda) cons
   vector<derivative> D;
   int hessianColsNbr = dynJacobianColsNbr*dynJacobianColsNbr;
   int thirdDerivativesColsNbr = hessianColsNbr*dynJacobianColsNbr;
-  for (third_derivatives_t::const_iterator it = third_derivatives.begin();
-       it != third_derivatives.end(); it++)
+  for (const auto & third_derivative : third_derivatives)
     {
-      int eq = it->first.first;
-      int var1 = it->first.second.first;
-      int var2 = it->first.second.second.first;
-      int var3 = it->first.second.second.second;
+      int eq = third_derivative.first.first;
+      int var1 = third_derivative.first.second.first;
+      int var2 = third_derivative.first.second.second.first;
+      int var3 = third_derivative.first.second.second.second;
 
       int id1 = getDynJacobianCol(var1);
       int id2 = getDynJacobianCol(var2);
@@ -6260,41 +6205,41 @@ DynamicModel::writeThirdDerivativesC_csr(const string &basename, bool cuda) cons
       vector<long unsigned int>  cols;
       long unsigned int col_nb = id1 * hessianColsNbr + id2 * dynJacobianColsNbr + id3;
       int thirdDColsNbr = hessianColsNbr*dynJacobianColsNbr;
-      derivative deriv(col_nb + eq *thirdDColsNbr, col_nb, eq, it->second);
+      derivative deriv(col_nb + eq *thirdDColsNbr, col_nb, eq, third_derivative.second);
       D.push_back(deriv);
       cols.push_back(col_nb);
       col_nb = id1 * hessianColsNbr + id3 * dynJacobianColsNbr + id2;
       if (find(cols.begin(), cols.end(), col_nb) == cols.end())
         {
-          derivative deriv(col_nb + eq *thirdDerivativesColsNbr, col_nb, eq, it->second);
+          derivative deriv(col_nb + eq *thirdDerivativesColsNbr, col_nb, eq, third_derivative.second);
           D.push_back(deriv);
           cols.push_back(col_nb);
         }
       col_nb = id2 * hessianColsNbr + id1 * dynJacobianColsNbr + id3;
       if (find(cols.begin(), cols.end(), col_nb) == cols.end())
         {
-          derivative deriv(col_nb + eq *thirdDerivativesColsNbr, col_nb, eq, it->second);
+          derivative deriv(col_nb + eq *thirdDerivativesColsNbr, col_nb, eq, third_derivative.second);
           D.push_back(deriv);
           cols.push_back(col_nb);
         }
       col_nb = id2 * hessianColsNbr + id3 * dynJacobianColsNbr + id1;
       if (find(cols.begin(), cols.end(), col_nb) == cols.end())
         {
-          derivative deriv(col_nb + eq *thirdDerivativesColsNbr, col_nb, eq, it->second);
+          derivative deriv(col_nb + eq *thirdDerivativesColsNbr, col_nb, eq, third_derivative.second);
           D.push_back(deriv);
           cols.push_back(col_nb);
         }
       col_nb = id3 * hessianColsNbr + id1 * dynJacobianColsNbr + id2;
       if (find(cols.begin(), cols.end(), col_nb) == cols.end())
         {
-          derivative deriv(col_nb + eq *thirdDerivativesColsNbr, col_nb, eq, it->second);
+          derivative deriv(col_nb + eq *thirdDerivativesColsNbr, col_nb, eq, third_derivative.second);
           D.push_back(deriv);
           cols.push_back(col_nb);
         }
       col_nb = id3 * hessianColsNbr + id2 * dynJacobianColsNbr + id1;
       if (find(cols.begin(), cols.end(), col_nb) == cols.end())
         {
-          derivative deriv(col_nb + eq *thirdDerivativesColsNbr, col_nb, eq, it->second);
+          derivative deriv(col_nb + eq *thirdDerivativesColsNbr, col_nb, eq, third_derivative.second);
           D.push_back(deriv);
         }
     }
@@ -6319,9 +6264,9 @@ DynamicModel::writeThirdDerivativesC_csr(const string &basename, bool cuda) cons
   // row_ptr must point to the relative address of the first element of the row
   int cumsum = 0;
   mDynamicModelFile << "row_ptr = [ 0";
-  for (vector<int>::iterator it = row_ptr.begin(); it != row_ptr.end(); ++it)
+  for (int & it : row_ptr)
     {
-      cumsum += *it;
+      cumsum += it;
       mDynamicModelFile << ", " << cumsum;
     }
   mDynamicModelFile << "];" << endl;
diff --git a/src/ExprNode.cc b/src/ExprNode.cc
index eadc2ae1..4ca81202 100644
--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
@@ -141,9 +141,8 @@ ExprNode::collectEndogenous(set<pair<int, int> > &result) const
 {
   set<pair<int, int> > symb_ids;
   collectDynamicVariables(eEndogenous, symb_ids);
-  for (set<pair<int, int> >::const_iterator it = symb_ids.begin();
-       it != symb_ids.end(); it++)
-    result.insert(make_pair(datatree.symbol_table.getTypeSpecificID(it->first), it->second));
+  for (const auto & symb_id : symb_ids)
+    result.insert(make_pair(datatree.symbol_table.getTypeSpecificID(symb_id.first), symb_id.second));
 }
 
 void
@@ -151,9 +150,8 @@ ExprNode::collectExogenous(set<pair<int, int> > &result) const
 {
   set<pair<int, int> > symb_ids;
   collectDynamicVariables(eExogenous, symb_ids);
-  for (set<pair<int, int> >::const_iterator it = symb_ids.begin();
-       it != symb_ids.end(); it++)
-    result.insert(make_pair(datatree.symbol_table.getTypeSpecificID(it->first), it->second));
+  for (const auto & symb_id : symb_ids)
+    result.insert(make_pair(datatree.symbol_table.getTypeSpecificID(symb_id.first), symb_id.second));
 }
 
 void
@@ -2015,9 +2013,8 @@ int
 UnaryOpNode::cost(const map<NodeTreeReference, temporary_terms_t> &temp_terms_map, bool is_matlab) const
 {
   // For a temporary term, the cost is null
-  for (map<NodeTreeReference, temporary_terms_t>::const_iterator it = temp_terms_map.begin();
-       it != temp_terms_map.end(); it++)
-    if (it->second.find(const_cast<UnaryOpNode *>(this)) != it->second.end())
+  for (const auto & it : temp_terms_map)
+    if (it.second.find(const_cast<UnaryOpNode *>(this)) != it.second.end())
       return 0;
 
   return cost(arg->cost(temp_terms_map, is_matlab), is_matlab);
@@ -3717,9 +3714,8 @@ int
 BinaryOpNode::cost(const map<NodeTreeReference, temporary_terms_t> &temp_terms_map, bool is_matlab) const
 {
   // For a temporary term, the cost is null
-  for (map<NodeTreeReference, temporary_terms_t>::const_iterator it = temp_terms_map.begin();
-       it != temp_terms_map.end(); it++)
-    if (it->second.find(const_cast<BinaryOpNode *>(this)) != it->second.end())
+  for (const auto & it : temp_terms_map)
+    if (it.second.find(const_cast<BinaryOpNode *>(this)) != it.second.end())
       return 0;
 
   int arg_cost = arg1->cost(temp_terms_map, is_matlab) + arg2->cost(temp_terms_map, is_matlab);
@@ -5254,9 +5250,8 @@ int
 TrinaryOpNode::cost(const map<NodeTreeReference, temporary_terms_t> &temp_terms_map, bool is_matlab) const
 {
   // For a temporary term, the cost is null
-  for (map<NodeTreeReference, temporary_terms_t>::const_iterator it = temp_terms_map.begin();
-       it != temp_terms_map.end(); it++)
-    if (it->second.find(const_cast<TrinaryOpNode *>(this)) != it->second.end())
+  for (const auto & it : temp_terms_map)
+    if (it.second.find(const_cast<TrinaryOpNode *>(this)) != it.second.end())
       return 0;
 
   int arg_cost = arg1->cost(temp_terms_map, is_matlab)
@@ -5986,8 +5981,8 @@ AbstractExternalFunctionNode::prepareForDerivation()
   if (preparedForDerivation)
     return;
 
-  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
-    (*it)->prepareForDerivation();
+  for (auto argument : arguments)
+    argument->prepareForDerivation();
 
   non_null_derivatives = arguments.at(0)->non_null_derivatives;
   for (int i = 1; i < (int) arguments.size(); i++)
@@ -6005,8 +6000,8 @@ AbstractExternalFunctionNode::computeDerivative(int deriv_id)
 {
   assert(datatree.external_functions_table.getNargs(symb_id) > 0);
   vector<expr_t> dargs;
-  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
-    dargs.push_back((*it)->getDerivative(deriv_id));
+  for (auto argument : arguments)
+    dargs.push_back(argument->getDerivative(deriv_id));
   return composeDerivatives(dargs);
 }
 
@@ -6015,9 +6010,8 @@ AbstractExternalFunctionNode::getChainRuleDerivative(int deriv_id, const map<int
 {
   assert(datatree.external_functions_table.getNargs(symb_id) > 0);
   vector<expr_t> dargs;
-  for (vector<expr_t>::const_iterator it = arguments.begin();
-       it != arguments.end(); it++)
-    dargs.push_back((*it)->getChainRuleDerivative(deriv_id, recursive_variables));
+  for (auto argument : arguments)
+    dargs.push_back(argument->getChainRuleDerivative(deriv_id, recursive_variables));
   return composeDerivatives(dargs);
 }
 
@@ -6027,9 +6021,8 @@ AbstractExternalFunctionNode::compileExternalFunctionArguments(ostream &CompileC
                                                                const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
                                                                const deriv_node_temp_terms_t &tef_terms) const
 {
-  for (vector<expr_t>::const_iterator it = arguments.begin();
-       it != arguments.end(); it++)
-    (*it)->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx,
+  for (auto argument : arguments)
+    argument->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx,
                    dynamic, steady_dynamic, tef_terms);
   return (arguments.size());
 }
@@ -6044,9 +6037,8 @@ AbstractExternalFunctionNode::collectVARLHSVariable(set<expr_t> &result) const
 void
 AbstractExternalFunctionNode::collectDynamicVariables(SymbolType type_arg, set<pair<int, int> > &result) const
 {
-  for (vector<expr_t>::const_iterator it = arguments.begin();
-       it != arguments.end(); it++)
-    (*it)->collectDynamicVariables(type_arg, result);
+  for (auto argument : arguments)
+    argument->collectDynamicVariables(type_arg, result);
 }
 
 void
@@ -6057,9 +6049,8 @@ AbstractExternalFunctionNode::collectTemporary_terms(const temporary_terms_t &te
     temporary_terms_inuse.insert(idx);
   else
     {
-      for (vector<expr_t>::const_iterator it = arguments.begin();
-           it != arguments.end(); it++)
-        (*it)->collectTemporary_terms(temporary_terms, temporary_terms_inuse, Curr_Block);
+      for (auto argument : arguments)
+        argument->collectTemporary_terms(temporary_terms, temporary_terms_inuse, Curr_Block);
     }
 }
 
@@ -6073,9 +6064,8 @@ int
 AbstractExternalFunctionNode::maxEndoLead() const
 {
   int val = 0;
-  for (vector<expr_t>::const_iterator it = arguments.begin();
-       it != arguments.end(); it++)
-    val = max(val, (*it)->maxEndoLead());
+  for (auto argument : arguments)
+    val = max(val, argument->maxEndoLead());
   return val;
 }
 
@@ -6083,9 +6073,8 @@ int
 AbstractExternalFunctionNode::maxExoLead() const
 {
   int val = 0;
-  for (vector<expr_t>::const_iterator it = arguments.begin();
-       it != arguments.end(); it++)
-    val = max(val, (*it)->maxExoLead());
+  for (auto argument : arguments)
+    val = max(val, argument->maxExoLead());
   return val;
 }
 
@@ -6093,9 +6082,8 @@ int
 AbstractExternalFunctionNode::maxEndoLag() const
 {
   int val = 0;
-  for (vector<expr_t>::const_iterator it = arguments.begin();
-       it != arguments.end(); it++)
-    val = max(val, (*it)->maxEndoLag());
+  for (auto argument : arguments)
+    val = max(val, argument->maxEndoLag());
   return val;
 }
 
@@ -6103,9 +6091,8 @@ int
 AbstractExternalFunctionNode::maxExoLag() const
 {
   int val = 0;
-  for (vector<expr_t>::const_iterator it = arguments.begin();
-       it != arguments.end(); it++)
-    val = max(val, (*it)->maxExoLag());
+  for (auto argument : arguments)
+    val = max(val, argument->maxExoLag());
   return val;
 }
 
@@ -6113,9 +6100,8 @@ int
 AbstractExternalFunctionNode::maxLead() const
 {
   int val = 0;
-  for (vector<expr_t>::const_iterator it = arguments.begin();
-       it != arguments.end(); it++)
-    val = max(val, (*it)->maxLead());
+  for (auto argument : arguments)
+    val = max(val, argument->maxLead());
   return val;
 }
 
@@ -6123,9 +6109,8 @@ int
 AbstractExternalFunctionNode::maxLag() const
 {
   int val = 0;
-  for (vector<expr_t>::const_iterator it = arguments.begin();
-       it != arguments.end(); it++)
-    val = max(val, (*it)->maxLag());
+  for (auto argument : arguments)
+    val = max(val, argument->maxLag());
   return val;
 }
 
@@ -6133,8 +6118,8 @@ expr_t
 AbstractExternalFunctionNode::undiff() const
 {
   vector<expr_t> arguments_subst;
-  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
-    arguments_subst.push_back((*it)->undiff());
+  for (auto argument : arguments)
+    arguments_subst.push_back(argument->undiff());
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
@@ -6142,27 +6127,24 @@ int
 AbstractExternalFunctionNode::VarMinLag() const
 {
 int val = 0;
-  for (vector<expr_t>::const_iterator it = arguments.begin();
-       it != arguments.end(); it++)
-    val = min(val, (*it)->VarMinLag());
+  for (auto argument : arguments)
+    val = min(val, argument->VarMinLag());
   return val;
 }
 
 void
 AbstractExternalFunctionNode::VarMaxLag(DataTree &static_datatree, set<expr_t> &static_lhs, int &max_lag) const
 {
-  for (vector<expr_t>::const_iterator it = arguments.begin();
-       it != arguments.end(); it++)
-    (*it)->VarMaxLag(static_datatree, static_lhs, max_lag);
+  for (auto argument : arguments)
+    argument->VarMaxLag(static_datatree, static_lhs, max_lag);
 }
 
 int
 AbstractExternalFunctionNode::PacMaxLag(vector<int> &lhs) const
 {
   int val = 0;
-  for (vector<expr_t>::const_iterator it = arguments.begin();
-       it != arguments.end(); it++)
-    val = max(val, (*it)->PacMaxLag(lhs));
+  for (auto argument : arguments)
+    val = max(val, argument->PacMaxLag(lhs));
   return val;
 }
 
@@ -6170,8 +6152,8 @@ expr_t
 AbstractExternalFunctionNode::decreaseLeadsLags(int n) const
 {
   vector<expr_t> arguments_subst;
-  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
-    arguments_subst.push_back((*it)->decreaseLeadsLags(n));
+  for (auto argument : arguments)
+    arguments_subst.push_back(argument->decreaseLeadsLags(n));
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
@@ -6179,8 +6161,8 @@ expr_t
 AbstractExternalFunctionNode::decreaseLeadsLagsPredeterminedVariables() const
 {
   vector<expr_t> arguments_subst;
-  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
-    arguments_subst.push_back((*it)->decreaseLeadsLagsPredeterminedVariables());
+  for (auto argument : arguments)
+    arguments_subst.push_back(argument->decreaseLeadsLagsPredeterminedVariables());
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
@@ -6188,8 +6170,8 @@ expr_t
 AbstractExternalFunctionNode::substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const
 {
   vector<expr_t> arguments_subst;
-  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
-    arguments_subst.push_back((*it)->substituteEndoLeadGreaterThanTwo(subst_table, neweqs, deterministic_model));
+  for (auto argument : arguments)
+    arguments_subst.push_back(argument->substituteEndoLeadGreaterThanTwo(subst_table, neweqs, deterministic_model));
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
@@ -6197,8 +6179,8 @@ expr_t
 AbstractExternalFunctionNode::substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
   vector<expr_t> arguments_subst;
-  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
-    arguments_subst.push_back((*it)->substituteEndoLagGreaterThanTwo(subst_table, neweqs));
+  for (auto argument : arguments)
+    arguments_subst.push_back(argument->substituteEndoLagGreaterThanTwo(subst_table, neweqs));
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
@@ -6206,8 +6188,8 @@ expr_t
 AbstractExternalFunctionNode::substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const
 {
   vector<expr_t> arguments_subst;
-  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
-    arguments_subst.push_back((*it)->substituteExoLead(subst_table, neweqs, deterministic_model));
+  for (auto argument : arguments)
+    arguments_subst.push_back(argument->substituteExoLead(subst_table, neweqs, deterministic_model));
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
@@ -6215,8 +6197,8 @@ expr_t
 AbstractExternalFunctionNode::substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
   vector<expr_t> arguments_subst;
-  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
-    arguments_subst.push_back((*it)->substituteExoLag(subst_table, neweqs));
+  for (auto argument : arguments)
+    arguments_subst.push_back(argument->substituteExoLag(subst_table, neweqs));
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
@@ -6224,8 +6206,8 @@ expr_t
 AbstractExternalFunctionNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const
 {
   vector<expr_t> arguments_subst;
-  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
-    arguments_subst.push_back((*it)->substituteExpectation(subst_table, neweqs, partial_information_model));
+  for (auto argument : arguments)
+    arguments_subst.push_back(argument->substituteExpectation(subst_table, neweqs, partial_information_model));
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
@@ -6233,23 +6215,23 @@ expr_t
 AbstractExternalFunctionNode::substituteAdl() const
 {
   vector<expr_t> arguments_subst;
-  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
-    arguments_subst.push_back((*it)->substituteAdl());
+  for (auto argument : arguments)
+    arguments_subst.push_back(argument->substituteAdl());
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
 void
 AbstractExternalFunctionNode::findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const
 {
-  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
-    (*it)->findDiffNodes(static_datatree, diff_table);
+  for (auto argument : arguments)
+    argument->findDiffNodes(static_datatree, diff_table);
 }
 
 void
 AbstractExternalFunctionNode::findUnaryOpNodesForAuxVarCreation(DataTree &static_datatree, diff_table_t &nodes) const
 {
-  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
-    (*it)->findUnaryOpNodesForAuxVarCreation(static_datatree, nodes);
+  for (auto argument : arguments)
+    argument->findUnaryOpNodesForAuxVarCreation(static_datatree, nodes);
 }
 
 expr_t
@@ -6257,8 +6239,8 @@ AbstractExternalFunctionNode::substituteDiff(DataTree &static_datatree, diff_tab
                                              vector<BinaryOpNode *> &neweqs) const
 {
   vector<expr_t> arguments_subst;
-  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
-    arguments_subst.push_back((*it)->substituteDiff(static_datatree, diff_table, subst_table, neweqs));
+  for (auto argument : arguments)
+    arguments_subst.push_back(argument->substituteDiff(static_datatree, diff_table, subst_table, neweqs));
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
@@ -6266,8 +6248,8 @@ expr_t
 AbstractExternalFunctionNode::substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
   vector<expr_t> arguments_subst;
-  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
-    arguments_subst.push_back((*it)->substituteUnaryOpNodes(static_datatree, nodes, subst_table, neweqs));
+  for (auto argument : arguments)
+    arguments_subst.push_back(argument->substituteUnaryOpNodes(static_datatree, nodes, subst_table, neweqs));
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
@@ -6275,8 +6257,8 @@ int
 AbstractExternalFunctionNode::countDiffs() const
 {
   int ndiffs = 0;
-  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
-    ndiffs += (*it)->countDiffs();
+  for (auto argument : arguments)
+    ndiffs += argument->countDiffs();
   return ndiffs;
 }
 
@@ -6284,8 +6266,8 @@ expr_t
 AbstractExternalFunctionNode::substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
 {
   vector<expr_t> arguments_subst;
-  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
-    arguments_subst.push_back((*it)->substitutePacExpectation(subst_table));
+  for (auto argument : arguments)
+    arguments_subst.push_back(argument->substitutePacExpectation(subst_table));
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
@@ -6293,8 +6275,8 @@ expr_t
 AbstractExternalFunctionNode::differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
   vector<expr_t> arguments_subst;
-  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
-    arguments_subst.push_back((*it)->differentiateForwardVars(subset, subst_table, neweqs));
+  for (auto argument : arguments)
+    arguments_subst.push_back(argument->differentiateForwardVars(subset, subst_table, neweqs));
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
@@ -6362,17 +6344,16 @@ bool
 AbstractExternalFunctionNode::containsEndogenous(void) const
 {
   bool result = false;
-  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
-    result = result || (*it)->containsEndogenous();
+  for (auto argument : arguments)
+    result = result || argument->containsEndogenous();
   return result;
 }
 
 bool
 AbstractExternalFunctionNode::containsExogenous() const
 {
-  for (vector<expr_t>::const_iterator it = arguments.begin();
-       it != arguments.end(); it++)
-    if ((*it)->containsExogenous())
+  for (auto argument : arguments)
+    if (argument->containsExogenous())
       return true;
   return false;
 }
@@ -6381,8 +6362,8 @@ expr_t
 AbstractExternalFunctionNode::replaceTrendVar() const
 {
   vector<expr_t> arguments_subst;
-  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
-    arguments_subst.push_back((*it)->replaceTrendVar());
+  for (auto argument : arguments)
+    arguments_subst.push_back(argument->replaceTrendVar());
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
@@ -6390,8 +6371,8 @@ expr_t
 AbstractExternalFunctionNode::detrend(int symb_id, bool log_trend, expr_t trend) const
 {
   vector<expr_t> arguments_subst;
-  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
-    arguments_subst.push_back((*it)->detrend(symb_id, log_trend, trend));
+  for (auto argument : arguments)
+    arguments_subst.push_back(argument->detrend(symb_id, log_trend, trend));
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
@@ -6399,16 +6380,16 @@ expr_t
 AbstractExternalFunctionNode::removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const
 {
   vector<expr_t> arguments_subst;
-  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
-    arguments_subst.push_back((*it)->removeTrendLeadLag(trend_symbols_map));
+  for (auto argument : arguments)
+    arguments_subst.push_back(argument->removeTrendLeadLag(trend_symbols_map));
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
 bool
 AbstractExternalFunctionNode::isInStaticForm() const
 {
-  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); ++it)
-    if (!(*it)->isInStaticForm())
+  for (auto argument : arguments)
+    if (!argument->isInStaticForm())
       return false;
   return true;
 }
@@ -6416,36 +6397,36 @@ AbstractExternalFunctionNode::isInStaticForm() const
 void
 AbstractExternalFunctionNode::setVarExpectationIndex(map<string, pair<SymbolList, int> > &var_model_info)
 {
-  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
-    (*it)->setVarExpectationIndex(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
 {
-  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
-    (*it)->walkPacParameters(pac_encountered, lhs, ec_params_and_vars, ar_params_and_vars);
+  for (auto argument : arguments)
+    argument->walkPacParameters(pac_encountered, lhs, ec_params_and_vars, ar_params_and_vars);
 }
 
 void
 AbstractExternalFunctionNode::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)
 {
-  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
-    (*it)->addParamInfoToPac(lhs_arg, ec_params_and_vars_arg, ar_params_and_vars_arg);
+  for (auto argument : arguments)
+    argument->addParamInfoToPac(lhs_arg, ec_params_and_vars_arg, ar_params_and_vars_arg);
 }
 
 void
 AbstractExternalFunctionNode::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)
 {
-  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
-    (*it)->fillPacExpectationVarInfo(model_name_arg, lhs_arg, max_lag_arg, nonstationary_arg, growth_symb_id_arg, equation_number_arg);
+  for (auto argument : arguments)
+    argument->fillPacExpectationVarInfo(model_name_arg, lhs_arg, max_lag_arg, nonstationary_arg, growth_symb_id_arg, equation_number_arg);
 }
 
 bool
 AbstractExternalFunctionNode::isVarModelReferenced(const string &model_info_name) const
 {
-  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); ++it)
-    if (!(*it)->isVarModelReferenced(model_info_name))
+  for (auto argument : arguments)
+    if (!argument->isVarModelReferenced(model_info_name))
       return true;
   return false;
 }
@@ -6453,8 +6434,8 @@ AbstractExternalFunctionNode::isVarModelReferenced(const string &model_info_name
 void
 AbstractExternalFunctionNode::getEndosAndMaxLags(map<string, int> &model_endos_and_lags) const
 {
-  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); ++it)
-    (*it)->getEndosAndMaxLags(model_endos_and_lags);
+  for (auto argument : arguments)
+    argument->getEndosAndMaxLags(model_endos_and_lags);
 }
 
 pair<int, expr_t>
@@ -6463,10 +6444,9 @@ AbstractExternalFunctionNode::normalizeEquation(int var_endo, vector<pair<int, p
   vector<pair<bool, expr_t> > V_arguments;
   vector<expr_t> V_expr_t;
   bool present = false;
-  for (vector<expr_t>::const_iterator it = arguments.begin();
-       it != arguments.end(); it++)
+  for (auto argument : arguments)
     {
-      V_arguments.push_back((*it)->normalizeEquation(var_endo, List_of_Op_RHS));
+      V_arguments.push_back(argument->normalizeEquation(var_endo, List_of_Op_RHS));
       present = present || V_arguments[V_arguments.size()-1].first;
       V_expr_t.push_back(V_arguments[V_arguments.size()-1].second);
     }
@@ -6516,11 +6496,10 @@ AbstractExternalFunctionNode::writePrhs(ostream &output, ExprNodeOutputType outp
 {
   output << "mxArray *prhs"<< ending << "[nrhs"<< ending << "];" << endl;
   int i = 0;
-  for (vector<expr_t>::const_iterator it = arguments.begin();
-       it != arguments.end(); it++)
+  for (auto argument : arguments)
     {
       output << "prhs" << ending << "[" << i++ << "] = mxCreateDoubleScalar("; // All external_function arguments are scalars
-      (*it)->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
+      argument->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
       output << ");" << endl;
     }
 }
@@ -6535,8 +6514,8 @@ expr_t
 AbstractExternalFunctionNode::substituteStaticAuxiliaryVariable() const
 {
   vector<expr_t> arguments_subst;
-  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
-    arguments_subst.push_back((*it)->substituteStaticAuxiliaryVariable());
+  for (auto argument : arguments)
+    arguments_subst.push_back(argument->substituteStaticAuxiliaryVariable());
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
@@ -6622,9 +6601,8 @@ ExternalFunctionNode::compileExternalFunctionOutput(ostream &CompileCode, unsign
   int first_deriv_symb_id = datatree.external_functions_table.getFirstDerivSymbID(symb_id);
   assert(first_deriv_symb_id != eExtFunSetButNoNameProvided);
 
-  for (vector<expr_t>::const_iterator it = arguments.begin();
-       it != arguments.end(); it++)
-    (*it)->compileExternalFunctionOutput(CompileCode, instruction_number, lhs_rhs, temporary_terms,
+  for (auto argument : arguments)
+    argument->compileExternalFunctionOutput(CompileCode, instruction_number, lhs_rhs, temporary_terms,
                                          map_idx, dynamic, steady_dynamic, tef_terms);
 
   if (!alreadyWrittenAsTefTerm(symb_id, tef_terms))
@@ -6717,9 +6695,8 @@ ExternalFunctionNode::writeExternalFunctionOutput(ostream &output, ExprNodeOutpu
   int first_deriv_symb_id = datatree.external_functions_table.getFirstDerivSymbID(symb_id);
   assert(first_deriv_symb_id != eExtFunSetButNoNameProvided);
 
-  for (vector<expr_t>::const_iterator it = arguments.begin();
-       it != arguments.end(); it++)
-    (*it)->writeExternalFunctionOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
+  for (auto argument : arguments)
+    argument->writeExternalFunctionOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
 
   if (!alreadyWrittenAsTefTerm(symb_id, tef_terms))
     {
@@ -6795,9 +6772,8 @@ ExternalFunctionNode::writeJsonExternalFunctionOutput(vector<string> &efout,
   int first_deriv_symb_id = datatree.external_functions_table.getFirstDerivSymbID(symb_id);
   assert(first_deriv_symb_id != eExtFunSetButNoNameProvided);
 
-  for (vector<expr_t>::const_iterator it = arguments.begin();
-       it != arguments.end(); it++)
-    (*it)->writeJsonExternalFunctionOutput(efout, temporary_terms, tef_terms, isdynamic);
+  for (auto argument : arguments)
+    argument->writeJsonExternalFunctionOutput(efout, temporary_terms, tef_terms, isdynamic);
 
   if (!alreadyWrittenAsTefTerm(symb_id, tef_terms))
     {
@@ -6827,9 +6803,8 @@ expr_t
 ExternalFunctionNode::toStatic(DataTree &static_datatree) const
 {
   vector<expr_t> static_arguments;
-  for (vector<expr_t>::const_iterator it = arguments.begin();
-       it != arguments.end(); it++)
-    static_arguments.push_back((*it)->toStatic(static_datatree));
+  for (auto argument : arguments)
+    static_arguments.push_back(argument->toStatic(static_datatree));
   return static_datatree.AddExternalFunction(symb_id, static_arguments);
 }
 
@@ -6837,18 +6812,16 @@ void
 ExternalFunctionNode::computeXrefs(EquationInfo &ei) const
 {
   vector<expr_t> dynamic_arguments;
-  for (vector<expr_t>::const_iterator it = arguments.begin();
-       it != arguments.end(); it++)
-    (*it)->computeXrefs(ei);
+  for (auto argument : arguments)
+    argument->computeXrefs(ei);
 }
 
 expr_t
 ExternalFunctionNode::cloneDynamic(DataTree &dynamic_datatree) const
 {
   vector<expr_t> dynamic_arguments;
-  for (vector<expr_t>::const_iterator it = arguments.begin();
-       it != arguments.end(); it++)
-    dynamic_arguments.push_back((*it)->cloneDynamic(dynamic_datatree));
+  for (auto argument : arguments)
+    dynamic_arguments.push_back(argument->cloneDynamic(dynamic_datatree));
   return dynamic_datatree.AddExternalFunction(symb_id, dynamic_arguments);
 }
 
@@ -7054,13 +7027,12 @@ FirstDerivExternalFunctionNode::writeExternalFunctionOutput(ostream &output, Exp
                << "prhs" << ending.str() << "[2] = mxCreateCellArray(2, dims" << ending.str() << ");"<< endl;
 
         int i = 0;
-        for (vector<expr_t>::const_iterator it = arguments.begin();
-             it != arguments.end(); it++)
+        for (auto argument : arguments)
           {
             output << "mxSetCell(prhs" << ending.str() << "[2], "
                    << i++ << ", "
                    << "mxCreateDoubleScalar("; // All external_function arguments are scalars
-            (*it)->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
+            argument->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
             output << "));" << endl;
           }
 
@@ -7206,9 +7178,8 @@ expr_t
 FirstDerivExternalFunctionNode::cloneDynamic(DataTree &dynamic_datatree) const
 {
   vector<expr_t> dynamic_arguments;
-  for (vector<expr_t>::const_iterator it = arguments.begin();
-       it != arguments.end(); it++)
-    dynamic_arguments.push_back((*it)->cloneDynamic(dynamic_datatree));
+  for (auto argument : arguments)
+    dynamic_arguments.push_back(argument->cloneDynamic(dynamic_datatree));
   return dynamic_datatree.AddFirstDerivExternalFunction(symb_id, dynamic_arguments,
                                                         inputIndex);
 }
@@ -7223,9 +7194,8 @@ expr_t
 FirstDerivExternalFunctionNode::toStatic(DataTree &static_datatree) const
 {
   vector<expr_t> static_arguments;
-  for (vector<expr_t>::const_iterator it = arguments.begin();
-       it != arguments.end(); it++)
-    static_arguments.push_back((*it)->toStatic(static_datatree));
+  for (auto argument : arguments)
+    static_arguments.push_back(argument->toStatic(static_datatree));
   return static_datatree.AddFirstDerivExternalFunction(symb_id, static_arguments,
                                                        inputIndex);
 }
@@ -7234,9 +7204,8 @@ void
 FirstDerivExternalFunctionNode::computeXrefs(EquationInfo &ei) const
 {
   vector<expr_t> dynamic_arguments;
-  for (vector<expr_t>::const_iterator it = arguments.begin();
-       it != arguments.end(); it++)
-    (*it)->computeXrefs(ei);
+  for (auto argument : arguments)
+    argument->computeXrefs(ei);
 }
 
 function<bool (expr_t)>
@@ -7415,13 +7384,12 @@ SecondDerivExternalFunctionNode::writeExternalFunctionOutput(ostream &output, Ex
                << "prhs" << ending.str() << "[3] = mxCreateCellArray(2, dims" << ending.str() << ");"<< endl;
 
         int i = 0;
-        for (vector<expr_t>::const_iterator it = arguments.begin();
-             it != arguments.end(); it++)
+        for (auto argument : arguments)
           {
             output << "mxSetCell(prhs" << ending.str() << "[3], "
                    << i++ << ", "
                    << "mxCreateDoubleScalar("; // All external_function arguments are scalars
-            (*it)->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
+            argument->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
             output << "));" << endl;
           }
 
@@ -7525,9 +7493,8 @@ expr_t
 SecondDerivExternalFunctionNode::cloneDynamic(DataTree &dynamic_datatree) const
 {
   vector<expr_t> dynamic_arguments;
-  for (vector<expr_t>::const_iterator it = arguments.begin();
-       it != arguments.end(); it++)
-    dynamic_arguments.push_back((*it)->cloneDynamic(dynamic_datatree));
+  for (auto argument : arguments)
+    dynamic_arguments.push_back(argument->cloneDynamic(dynamic_datatree));
   return dynamic_datatree.AddSecondDerivExternalFunction(symb_id, dynamic_arguments,
                                                          inputIndex1, inputIndex2);
 }
@@ -7542,9 +7509,8 @@ expr_t
 SecondDerivExternalFunctionNode::toStatic(DataTree &static_datatree) const
 {
   vector<expr_t> static_arguments;
-  for (vector<expr_t>::const_iterator it = arguments.begin();
-       it != arguments.end(); it++)
-    static_arguments.push_back((*it)->toStatic(static_datatree));
+  for (auto argument : arguments)
+    static_arguments.push_back(argument->toStatic(static_datatree));
   return static_datatree.AddSecondDerivExternalFunction(symb_id, static_arguments,
                                                         inputIndex1, inputIndex2);
 }
@@ -7553,9 +7519,8 @@ void
 SecondDerivExternalFunctionNode::computeXrefs(EquationInfo &ei) const
 {
   vector<expr_t> dynamic_arguments;
-  for (vector<expr_t>::const_iterator it = arguments.begin();
-       it != arguments.end(); it++)
-    (*it)->computeXrefs(ei);
+  for (auto argument : arguments)
+    argument->computeXrefs(ei);
 }
 
 void
diff --git a/src/ExternalFunctionsTable.hh b/src/ExternalFunctionsTable.hh
index 42d18391..6c2e6b5d 100644
--- a/src/ExternalFunctionsTable.hh
+++ b/src/ExternalFunctionsTable.hh
@@ -116,9 +116,8 @@ inline int
 ExternalFunctionsTable::get_total_number_of_unique_model_block_external_functions() const
 {
   int number_of_unique_model_block_external_functions = 0;
-  for (external_function_table_type::const_iterator it = externalFunctionTable.begin();
-       it != externalFunctionTable.end(); it++)
-    if (it->second.nargs > 0)
+  for (auto it : externalFunctionTable)
+    if (it.second.nargs > 0)
       number_of_unique_model_block_external_functions++;
 
   return number_of_unique_model_block_external_functions;
diff --git a/src/ModFile.cc b/src/ModFile.cc
index 873127df..1d3c52bd 100644
--- a/src/ModFile.cc
+++ b/src/ModFile.cc
@@ -48,9 +48,8 @@ ModFile::ModFile(WarningConsolidation &warnings_arg)
 
 ModFile::~ModFile()
 {
-  for (vector<Statement *>::iterator it = statements.begin();
-       it != statements.end(); it++)
-    delete (*it);
+  for (auto & statement : statements)
+    delete statement;
 }
 
 void
@@ -112,9 +111,8 @@ ModFile::addStatementAtFront(Statement *st)
 void
 ModFile::checkPass(bool nostrict, bool stochastic)
 {
-  for (vector<Statement *>::iterator it = statements.begin();
-       it != statements.end(); it++)
-    (*it)->checkPass(mod_file_struct, warnings);
+  for (auto & statement : statements)
+    statement->checkPass(mod_file_struct, warnings);
 
   // Check the steady state block
   steady_state_model.checkPass(mod_file_struct, warnings);
@@ -332,8 +330,8 @@ ModFile::checkPass(bool nostrict, bool stochastic)
   if (unusedExo.size() > 0)
     {
       ostringstream unused_exos;
-      for (set<int>::iterator it = unusedExo.begin(); it != unusedExo.end(); it++)
-        unused_exos << symbol_table.getName(*it) << " ";
+      for (int it : unusedExo)
+        unused_exos << symbol_table.getName(it) << " ";
 
       if (nostrict)
         warnings << "WARNING: " << unused_exos.str()
@@ -358,10 +356,10 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const
   if (nostrict)
     {
       set<int> unusedEndogs = dynamic_model.findUnusedEndogenous();
-      for (set<int>::iterator it = unusedEndogs.begin(); it != unusedEndogs.end(); it++)
+      for (int unusedEndog : unusedEndogs)
         {
-          symbol_table.changeType(*it, eUnusedEndogenous);
-          warnings << "WARNING: '" << symbol_table.getName(*it)
+          symbol_table.changeType(unusedEndog, eUnusedEndogenous);
+          warnings << "WARNING: '" << symbol_table.getName(unusedEndog)
                    << "' not used in model block, removed by nostrict command-line option" << endl;
         }
     }
@@ -462,9 +460,9 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const
   if (mod_file_struct.ramsey_model_present)
     {
       StaticModel *planner_objective = NULL;
-      for (vector<Statement *>::iterator it = statements.begin(); it != statements.end(); it++)
+      for (auto & statement : statements)
         {
-          PlannerObjectiveStatement *pos = dynamic_cast<PlannerObjectiveStatement *>(*it);
+          PlannerObjectiveStatement *pos = dynamic_cast<PlannerObjectiveStatement *>(statement);
           if (pos != NULL)
             planner_objective = pos->getPlannerObjective();
         }
@@ -597,9 +595,9 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const
     }
 
   if (mod_file_struct.ramsey_policy_present)
-    for (vector<Statement *>::iterator it = statements.begin(); it != statements.end(); it++)
+    for (auto & statement : statements)
       {
-        RamseyPolicyStatement *rps = dynamic_cast<RamseyPolicyStatement *>(*it);
+        RamseyPolicyStatement *rps = dynamic_cast<RamseyPolicyStatement *>(statement);
         if (rps != NULL)
           rps->checkRamseyPolicyList();
       }
@@ -727,9 +725,8 @@ ModFile::computingPass(bool no_tmp_terms, FileOutputType output, int params_deri
         }
     }
 
-  for (vector<Statement *>::iterator it = statements.begin();
-       it != statements.end(); it++)
-    (*it)->computingPass();
+  for (auto & statement : statements)
+    statement->computingPass();
 }
 
 void
@@ -848,14 +845,14 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool clear_glo
   if (parallel_local_files.size() > 0)
     {
       mOutputFile << "options_.parallel_info.local_files = {" << endl;
-      for (size_t i = 0; i < parallel_local_files.size(); i++)
+      for (const auto & parallel_local_file : parallel_local_files)
         {
-          size_t j = parallel_local_files[i].find_last_of("/\\");
+          size_t j = parallel_local_file.find_last_of("/\\");
           if (j == string::npos)
-            mOutputFile << "'', '" << parallel_local_files[i] << "';" << endl;
+            mOutputFile << "'', '" << parallel_local_file << "';" << endl;
           else
-            mOutputFile << "'" << parallel_local_files[i].substr(0, j+1) << "', '"
-                        << parallel_local_files[i].substr(j+1, string::npos) << "';" << endl;
+            mOutputFile << "'" << parallel_local_file.substr(0, j+1) << "', '"
+                        << parallel_local_file.substr(j+1, string::npos) << "';" << endl;
         }
       mOutputFile << "};" << endl;
     }
@@ -985,14 +982,13 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool clear_glo
   //      (*it)->writeOutput(mOutputFile, basename, minimal_workspace);
   //  dynamic_model.writeVarExpectationFunctions(mOutputFile);
 
-  for (vector<Statement *>::const_iterator it = statements.begin();
-       it != statements.end(); it++)
+  for (auto statement : statements)
     {
-      (*it)->writeOutput(mOutputFile, basename, minimal_workspace);
+      statement->writeOutput(mOutputFile, basename, minimal_workspace);
 
       /* Special treatment for initval block: insert initial values for the
          auxiliary variables and initialize exo det */
-      InitValStatement *ivs = dynamic_cast<InitValStatement *>(*it);
+      InitValStatement *ivs = dynamic_cast<InitValStatement *>(statement);
       if (ivs != NULL)
         {
           static_model.writeAuxVarInitval(mOutputFile, oMatlabOutsideModel);
@@ -1000,17 +996,17 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool clear_glo
         }
 
       // Special treatment for endval block: insert initial values for the auxiliary variables
-      EndValStatement *evs = dynamic_cast<EndValStatement *>(*it);
+      EndValStatement *evs = dynamic_cast<EndValStatement *>(statement);
       if (evs != NULL)
         static_model.writeAuxVarInitval(mOutputFile, oMatlabOutsideModel);
 
       // Special treatment for load params and steady state statement: insert initial values for the auxiliary variables
-      LoadParamsAndSteadyStateStatement *lpass = dynamic_cast<LoadParamsAndSteadyStateStatement *>(*it);
+      LoadParamsAndSteadyStateStatement *lpass = dynamic_cast<LoadParamsAndSteadyStateStatement *>(statement);
       if (lpass && !no_static)
         static_model.writeAuxVarInitval(mOutputFile, oMatlabOutsideModel);
 
       // Special treatement for Var Models
-      VarModelStatement *vms = dynamic_cast<VarModelStatement *>(*it);
+      VarModelStatement *vms = dynamic_cast<VarModelStatement *>(statement);
       if (vms != NULL)
         vms->createVarModelMFunction(mOutputFile, dynamic_model.getVarExpectationFunctionsToWrite());
     }
@@ -1162,9 +1158,8 @@ ModFile::writeModelC(const string &basename) const
     }
 
   // Print statements
-  for (vector<Statement *>::const_iterator it = statements.begin();
-       it != statements.end(); it++)
-    (*it)->writeCOutput(mDriverCFile, basename);
+  for (auto statement : statements)
+    statement->writeCOutput(mDriverCFile, basename);
 
   mDriverCFile << "} DynareInfo;" << endl;
   mDriverCFile.close();
@@ -1267,9 +1262,8 @@ ModFile::writeModelCC(const string &basename) const
     }
 
   // Print statements
-  for (vector<Statement *>::const_iterator it = statements.begin();
-       it != statements.end(); it++)
-    (*it)->writeCOutput(mDriverCFile, basename);
+  for (auto statement : statements)
+    statement->writeCOutput(mDriverCFile, basename);
 
   mDriverCFile << "};" << endl;
   mDriverCFile.close();
@@ -1404,9 +1398,8 @@ ModFile::writeExternalFilesJulia(const string &basename, FileOutputType output,
   steady_state_model.writeSteadyStateFile(basename, mod_file_struct.ramsey_model_present, true);
 
   // Print statements (includes parameter values)
-  for (vector<Statement *>::const_iterator it = statements.begin();
-       it != statements.end(); it++)
-    (*it)->writeJuliaOutput(jlOutputFile, basename);
+  for (auto statement : statements)
+    statement->writeJuliaOutput(jlOutputFile, basename);
 
   jlOutputFile << "model_.static = " << basename << "Static.static!" << endl
                << "model_.dynamic = " << basename << "Dynamic.dynamic!" << endl
diff --git a/src/ModelTree.cc b/src/ModelTree.cc
index 11c12a04..343c93e4 100644
--- a/src/ModelTree.cc
+++ b/src/ModelTree.cc
@@ -51,8 +51,8 @@ ModelTree::computeNormalization(const jacob_map_t &contemporaneous_jacobian, boo
   // Fill in the graph
   set<pair<int, int> > endo;
 
-  for (jacob_map_t::const_iterator it = contemporaneous_jacobian.begin(); it != contemporaneous_jacobian.end(); it++)
-    add_edge(it->first.first + n, it->first.second, g);
+  for (const auto & it : contemporaneous_jacobian)
+    add_edge(it.first.first + n, it.first.second, g);
 
   // Compute maximum cardinality matching
   vector<int> mate_map(2*n);
@@ -153,8 +153,8 @@ ModelTree::computeNonSingularNormalization(jacob_map_t &contemporaneous_jacobian
     if (fabs(iter->second) > max_val[iter->first.first])
       max_val[iter->first.first] = fabs(iter->second);
 
-  for (jacob_map_t::iterator iter = normalized_contemporaneous_jacobian.begin(); iter != normalized_contemporaneous_jacobian.end(); iter++)
-    iter->second /= max_val[iter->first.first];
+  for (auto & iter : normalized_contemporaneous_jacobian)
+    iter.second /= max_val[iter.first.first];
 
   //We start with the highest value of the cutoff and try to normalize the model
   double current_cutoff = 0.99999999;
@@ -164,9 +164,9 @@ ModelTree::computeNonSingularNormalization(jacob_map_t &contemporaneous_jacobian
     {
       jacob_map_t tmp_normalized_contemporaneous_jacobian;
       int suppress = 0;
-      for (jacob_map_t::iterator iter = normalized_contemporaneous_jacobian.begin(); iter != normalized_contemporaneous_jacobian.end(); iter++)
-        if (fabs(iter->second) > max(current_cutoff, cutoff))
-          tmp_normalized_contemporaneous_jacobian[make_pair(iter->first.first, iter->first.second)] = iter->second;
+      for (auto & iter : normalized_contemporaneous_jacobian)
+        if (fabs(iter.second) > max(current_cutoff, cutoff))
+          tmp_normalized_contemporaneous_jacobian[make_pair(iter.first.first, iter.first.second)] = iter.second;
         else
           suppress++;
 
@@ -192,8 +192,8 @@ ModelTree::computeNonSingularNormalization(jacob_map_t &contemporaneous_jacobian
         {
           endo.clear();
           equations[i]->collectEndogenous(endo);
-          for (set<pair<int, int> >::const_iterator it = endo.begin(); it != endo.end(); it++)
-            tmp_normalized_contemporaneous_jacobian[make_pair(i, it->first)] = 1;
+          for (const auto & it : endo)
+            tmp_normalized_contemporaneous_jacobian[make_pair(i, it.first)] = 1;
         }
       check = computeNormalization(tmp_normalized_contemporaneous_jacobian, true);
       if (check)
@@ -307,8 +307,8 @@ ModelTree::evaluateAndReduceJacobian(const eval_context_t &eval_context, jacob_m
     }
 
   // Get rid of the elements of the Jacobian matrix below the cutoff
-  for (set<pair<int, int> >::const_iterator it = jacobian_elements_to_delete.begin(); it != jacobian_elements_to_delete.end(); it++)
-    first_derivatives.erase(*it);
+  for (const auto & it : jacobian_elements_to_delete)
+    first_derivatives.erase(it);
 
   if (jacobian_elements_to_delete.size() > 0)
     {
@@ -340,13 +340,13 @@ ModelTree::computePrologueAndEpilogue(const jacob_map_t &static_jacobian_arg, ve
         {
           endo.clear();
           equations[i]->collectEndogenous(endo);
-          for (set<pair<int, int> >::const_iterator it = endo.begin(); it != endo.end(); it++)
-            IM[i * n + endo2eq[it->first]] = true;
+          for (const auto & it : endo)
+            IM[i * n + endo2eq[it.first]] = true;
         }
     }
   else
-    for (jacob_map_t::const_iterator it = static_jacobian_arg.begin(); it != static_jacobian_arg.end(); it++)
-      IM[it->first.first * n + endo2eq[it->first.second]] = true;
+    for (const auto & it : static_jacobian_arg)
+      IM[it.first.first * n + endo2eq[it.first.second]] = true;
   bool something_has_been_done = true;
   prologue = 0;
   int k = 0;
@@ -508,11 +508,11 @@ ModelTree::getVariableLeadLagByBlock(const dynamic_jacob_map_t &dynamic_jacobian
           equation_blck[equation_reordered[i]] = i- (nb_endo - nb_blck_sim - prologue - epilogue);
         }
     }
-  for (dynamic_jacob_map_t::const_iterator it = dynamic_jacobian.begin(); it != dynamic_jacobian.end(); it++)
+  for (const auto & it : dynamic_jacobian)
     {
-      int lag = it->first.first;
-      int j_1 = it->first.second.first;
-      int i_1 = it->first.second.second;
+      int lag = it.first.first;
+      int j_1 = it.first.second.first;
+      int i_1 = it.first.second.second;
       if (variable_blck[i_1] == equation_blck[j_1])
         {
           if (lag > variable_lead_lag[i_1].second)
@@ -555,8 +555,8 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
         {
           endo.clear();
           equations[i]->collectEndogenous(endo);
-          for (set<pair<int, int> >::const_iterator it = endo.begin(); it != endo.end(); it++)
-            tmp_normalized_contemporaneous_jacobian[make_pair(i, it->first)] = 1;
+          for (const auto & it : endo)
+            tmp_normalized_contemporaneous_jacobian[make_pair(i, it.first)] = 1;
 
         }
     }
@@ -674,33 +674,33 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
       //First we have the recursive equations conditional on feedback variables
       for (int j = 0; j < 4; j++)
         {
-          for (vector<int>::iterator its = Reordered_Vertice.begin(); its != Reordered_Vertice.end(); its++)
+          for (int & its : Reordered_Vertice)
             {
               bool something_done = false;
-              if      (j == 2 && variable_lag_lead[tmp_variable_reordered[*its +prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[*its +prologue]].second != 0)
+              if      (j == 2 && variable_lag_lead[tmp_variable_reordered[its +prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[its +prologue]].second != 0)
                 {
                   n_mixed[prologue+i]++;
                   something_done = true;
                 }
-              else if (j == 3 && variable_lag_lead[tmp_variable_reordered[*its +prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[*its +prologue]].second != 0)
+              else if (j == 3 && variable_lag_lead[tmp_variable_reordered[its +prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[its +prologue]].second != 0)
                 {
                   n_forward[prologue+i]++;
                   something_done = true;
                 }
-              else if (j == 1 && variable_lag_lead[tmp_variable_reordered[*its +prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[*its +prologue]].second == 0)
+              else if (j == 1 && variable_lag_lead[tmp_variable_reordered[its +prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[its +prologue]].second == 0)
                 {
                   n_backward[prologue+i]++;
                   something_done = true;
                 }
-              else if (j == 0 && variable_lag_lead[tmp_variable_reordered[*its +prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[*its +prologue]].second == 0)
+              else if (j == 0 && variable_lag_lead[tmp_variable_reordered[its +prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[its +prologue]].second == 0)
                 {
                   n_static[prologue+i]++;
                   something_done = true;
                 }
               if (something_done)
                 {
-                  equation_reordered[order] = tmp_equation_reordered[*its+prologue];
-                  variable_reordered[order] = tmp_variable_reordered[*its+prologue];
+                  equation_reordered[order] = tmp_equation_reordered[its+prologue];
+                  variable_reordered[order] = tmp_variable_reordered[its+prologue];
                   order++;
                 }
             }
@@ -709,33 +709,33 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
       //Second we have the equations related to the feedback variables
       for (int j = 0; j < 4; j++)
         {
-          for (set<int>::iterator its = feed_back_vertices.begin(); its != feed_back_vertices.end(); its++)
+          for (int feed_back_vertice : feed_back_vertices)
             {
               bool something_done = false;
-              if      (j == 2 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(*its, G)]+prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(*its, G)]+prologue]].second != 0)
+              if      (j == 2 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]].second != 0)
                 {
                   n_mixed[prologue+i]++;
                   something_done = true;
                 }
-              else if (j == 3 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(*its, G)]+prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(*its, G)]+prologue]].second != 0)
+              else if (j == 3 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]].second != 0)
                 {
                   n_forward[prologue+i]++;
                   something_done = true;
                 }
-              else if (j == 1 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(*its, G)]+prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(*its, G)]+prologue]].second == 0)
+              else if (j == 1 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]].second == 0)
                 {
                   n_backward[prologue+i]++;
                   something_done = true;
                 }
-              else if (j == 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(*its, G)]+prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(*its, G)]+prologue]].second == 0)
+              else if (j == 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]].second == 0)
                 {
                   n_static[prologue+i]++;
                   something_done = true;
                 }
               if (something_done)
                 {
-                  equation_reordered[order] = tmp_equation_reordered[v_index[vertex(*its, G)]+prologue];
-                  variable_reordered[order] = tmp_variable_reordered[v_index[vertex(*its, G)]+prologue];
+                  equation_reordered[order] = tmp_equation_reordered[v_index[vertex(feed_back_vertice, G)]+prologue];
+                  variable_reordered[order] = tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue];
                   order++;
                 }
             }
@@ -832,10 +832,10 @@ ModelTree::reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_j
         {
           endo.clear();
           equations[equation_reordered[count_equ]]->collectEndogenous(endo);
-          for (set<pair<int, int> >::const_iterator it = endo.begin(); it != endo.end(); it++)
+          for (const auto & it : endo)
             {
-              int curr_variable = it->first;
-              int curr_lag = it->second;
+              int curr_variable = it.first;
+              int curr_lag = it.second;
               vector<int>::const_iterator it1 = find(variable_reordered.begin()+first_count_equ, variable_reordered.begin()+(first_count_equ+Blck_Size), curr_variable);
               if (it1 != variable_reordered.begin()+(first_count_equ+Blck_Size))
                 if (dynamic_jacobian.find(make_pair(curr_lag, make_pair(equation_reordered[count_equ], curr_variable))) != dynamic_jacobian.end())
@@ -1010,8 +1010,8 @@ ModelTree::ModelTree(SymbolTable &symbol_table_arg,
   mfs(0)
 
 {
-  for (int i = 0; i < 3; i++)
-    NNZDerivatives[i] = 0;
+  for (int & NNZDerivative : NNZDerivatives)
+    NNZDerivative = 0;
 }
 
 int
@@ -1035,15 +1035,14 @@ ModelTree::writeDerivative(ostream &output, int eq, int symb_id, int lag,
 void
 ModelTree::computeJacobian(const set<int> &vars)
 {
-  for (set<int>::const_iterator it = vars.begin();
-       it != vars.end(); it++)
+  for (int var : vars)
     {
       for (int eq = 0; eq < (int) equations.size(); eq++)
         {
-          expr_t d1 = equations[eq]->getDerivative(*it);
+          expr_t d1 = equations[eq]->getDerivative(var);
           if (d1 == Zero)
             continue;
-          first_derivatives[make_pair(eq, *it)] = d1;
+          first_derivatives[make_pair(eq, var)] = d1;
           ++NNZDerivatives[0];
         }
     }
@@ -1060,10 +1059,8 @@ ModelTree::computeHessian(const set<int> &vars)
       expr_t d1 = it->second;
 
       // Store only second derivatives with var2 <= var1
-      for (set<int>::const_iterator it2 = vars.begin();
-           it2 != vars.end(); it2++)
+      for (int var2 : vars)
         {
-          int var2 = *it2;
           if (var2 > var1)
             continue;
 
@@ -1094,10 +1091,8 @@ ModelTree::computeThirdDerivatives(const set<int> &vars)
       expr_t d2 = it->second;
 
       // Store only third derivatives such that var3 <= var2 <= var1
-      for (set<int>::const_iterator it2 = vars.begin();
-           it2 != vars.end(); it2++)
+      for (int var3 : vars)
         {
-          int var3 = *it2;
           if (var3 > var2)
             continue;
 
@@ -1129,14 +1124,13 @@ ModelTree::computeTemporaryTerms(bool is_matlab)
   // Collect all model local variables appearing in equations. See #101
   // All used model local variables are automatically set as temporary variables
   set<int> used_local_vars;
-  for (size_t i = 0; i < equations.size(); i++)
-    equations[i]->collectVariables(eModelLocalVariable, used_local_vars);
+  for (auto & equation : equations)
+    equation->collectVariables(eModelLocalVariable, used_local_vars);
 
-  for (set<int>::const_iterator it = used_local_vars.begin();
-       it != used_local_vars.end(); it++)
+  for (int used_local_var : used_local_vars)
     {
-      VariableNode *v = AddVariable(*it);
-      temporary_terms_mlv[v] = local_variables_table.find(*it)->second;
+      VariableNode *v = AddVariable(used_local_var);
+      temporary_terms_mlv[v] = local_variables_table.find(used_local_var)->second;
       reference_count[v] = make_pair(MIN_COST(is_matlab)+1, eResiduals);
     }
 
@@ -1146,27 +1140,23 @@ ModelTree::computeTemporaryTerms(bool is_matlab)
   temp_terms_map[eSecondDeriv] = temporary_terms_g2;
   temp_terms_map[eThirdDeriv] = temporary_terms_g3;
 
-  for (vector<BinaryOpNode *>::iterator it = equations.begin();
-       it != equations.end(); it++)
-    (*it)->computeTemporaryTerms(reference_count,
+  for (auto & equation : equations)
+    equation->computeTemporaryTerms(reference_count,
                                  temp_terms_map,
                                  is_matlab, eResiduals);
 
-  for (first_derivatives_t::iterator it = first_derivatives.begin();
-       it != first_derivatives.end(); it++)
-    it->second->computeTemporaryTerms(reference_count,
+  for (auto & first_derivative : first_derivatives)
+    first_derivative.second->computeTemporaryTerms(reference_count,
                                       temp_terms_map,
                                       is_matlab, eFirstDeriv);
 
-  for (second_derivatives_t::iterator it = second_derivatives.begin();
-       it != second_derivatives.end(); it++)
-    it->second->computeTemporaryTerms(reference_count,
+  for (auto & second_derivative : second_derivatives)
+    second_derivative.second->computeTemporaryTerms(reference_count,
                                       temp_terms_map,
                                       is_matlab, eSecondDeriv);
 
-  for (third_derivatives_t::iterator it = third_derivatives.begin();
-       it != third_derivatives.end(); it++)
-    it->second->computeTemporaryTerms(reference_count,
+  for (auto & third_derivative : third_derivatives)
+    third_derivative.second->computeTemporaryTerms(reference_count,
                                       temp_terms_map,
                                       is_matlab, eThirdDeriv);
 
@@ -1184,21 +1174,17 @@ ModelTree::computeTemporaryTerms(bool is_matlab)
        it != temporary_terms_mlv.end(); it++)
     temporary_terms_idxs[it->first] = idx++;
 
-  for (temporary_terms_t::const_iterator it = temporary_terms_res.begin();
-       it != temporary_terms_res.end(); it++)
-    temporary_terms_idxs[*it] = idx++;
+  for (auto temporary_terms_re : temporary_terms_res)
+    temporary_terms_idxs[temporary_terms_re] = idx++;
 
-  for (temporary_terms_t::const_iterator it = temporary_terms_g1.begin();
-       it != temporary_terms_g1.end(); it++)
-    temporary_terms_idxs[*it] = idx++;
+  for (auto it : temporary_terms_g1)
+    temporary_terms_idxs[it] = idx++;
 
-  for (temporary_terms_t::const_iterator it = temporary_terms_g2.begin();
-       it != temporary_terms_g2.end(); it++)
-    temporary_terms_idxs[*it] = idx++;
+  for (auto it : temporary_terms_g2)
+    temporary_terms_idxs[it] = idx++;
 
-  for (temporary_terms_t::const_iterator it = temporary_terms_g3.begin();
-       it != temporary_terms_g3.end(); it++)
-    temporary_terms_idxs[*it] = idx++;
+  for (auto it : temporary_terms_g3)
+    temporary_terms_idxs[it] = idx++;
 }
 
 void
@@ -1207,24 +1193,23 @@ ModelTree::writeModelLocalVariableTemporaryTerms(const temporary_terms_t &tto, c
                                                  deriv_node_temp_terms_t &tef_terms) const
 {
   temporary_terms_t tt2;
-  for (map<expr_t, expr_t>::const_iterator it = tt.begin();
-       it != tt.end(); it++)
+  for (auto it : tt)
     {
       if (IS_C(output_type))
         output << "double ";
       else if (IS_JULIA(output_type))
         output << "    @inbounds const ";
 
-      it->first->writeOutput(output, output_type, tto, temporary_terms_idxs, tef_terms);
+      it.first->writeOutput(output, output_type, tto, temporary_terms_idxs, tef_terms);
       output << " = ";
-      it->second->writeOutput(output, output_type, tt2, temporary_terms_idxs, tef_terms);
+      it.second->writeOutput(output, output_type, tt2, temporary_terms_idxs, tef_terms);
 
       if (IS_C(output_type) || IS_MATLAB(output_type))
         output << ";";
       output << endl;
 
       // Insert current node into tt2
-      tt2.insert(it->first);
+      tt2.insert(it.first);
     }
 }
 
@@ -1269,16 +1254,15 @@ ModelTree::writeJsonTemporaryTerms(const temporary_terms_t &tt, const temporary_
   temporary_terms_t tt2 = ttm1;
 
   output << "\"external_functions_temporary_terms_" << concat << "\": [";
-  for (temporary_terms_t::const_iterator it = tt.begin();
-       it != tt.end(); it++)
-    if (ttm1.find(*it) == ttm1.end())
+  for (auto it : tt)
+    if (ttm1.find(it) == ttm1.end())
       {
-        if (dynamic_cast<AbstractExternalFunctionNode *>(*it) != NULL)
+        if (dynamic_cast<AbstractExternalFunctionNode *>(it) != NULL)
           {
             if (wrote_term)
               output << ", ";
             vector<string> efout;
-            (*it)->writeJsonExternalFunctionOutput(efout, tt2, tef_terms);
+            it->writeJsonExternalFunctionOutput(efout, tt2, tef_terms);
             for (vector<string>::const_iterator it1 = efout.begin(); it1 != efout.end(); it1++)
               {
                 if (it1 != efout.begin())
@@ -1287,7 +1271,7 @@ ModelTree::writeJsonTemporaryTerms(const temporary_terms_t &tt, const temporary_
               }
             wrote_term = true;
           }
-        tt2.insert(*it);
+        tt2.insert(it);
       }
 
   tt2 = ttm1;
@@ -1424,11 +1408,11 @@ bool
 ModelTree::testNestedParenthesis(const string &str) const
 {
   int open = 0;
-  for (size_t i = 0; i < str.length(); i++)
+  for (char i : str)
     {
-      if (str.at(i) == '(')
+      if (i == '(')
         open++;
-      else if (str.at(i) == ')')
+      else if (i == ')')
         open--;
       if (open > 32)
         return true;
@@ -1443,29 +1427,28 @@ ModelTree::compileTemporaryTerms(ostream &code_file, unsigned int &instruction_n
   temporary_terms_t tt2;
   // To store the functions that have already been written in the form TEF* = ext_fun();
   deriv_node_temp_terms_t tef_terms;
-  for (temporary_terms_t::const_iterator it = tt.begin();
-       it != tt.end(); it++)
+  for (auto it : tt)
     {
-      if (dynamic_cast<AbstractExternalFunctionNode *>(*it) != NULL)
+      if (dynamic_cast<AbstractExternalFunctionNode *>(it) != NULL)
         {
-          (*it)->compileExternalFunctionOutput(code_file, instruction_number, false, tt2, map_idx, dynamic, steady_dynamic, tef_terms);
+          it->compileExternalFunctionOutput(code_file, instruction_number, false, tt2, map_idx, dynamic, steady_dynamic, tef_terms);
         }
 
-      FNUMEXPR_ fnumexpr(TemporaryTerm, (int)(map_idx.find((*it)->idx)->second));
+      FNUMEXPR_ fnumexpr(TemporaryTerm, (int)(map_idx.find(it->idx)->second));
       fnumexpr.write(code_file, instruction_number);
-      (*it)->compile(code_file, instruction_number, false, tt2, map_idx, dynamic, steady_dynamic, tef_terms);
+      it->compile(code_file, instruction_number, false, tt2, map_idx, dynamic, steady_dynamic, tef_terms);
       if (dynamic)
         {
-          FSTPT_ fstpt((int)(map_idx.find((*it)->idx)->second));
+          FSTPT_ fstpt((int)(map_idx.find(it->idx)->second));
           fstpt.write(code_file, instruction_number);
         }
       else
         {
-          FSTPST_ fstpst((int)(map_idx.find((*it)->idx)->second));
+          FSTPST_ fstpst((int)(map_idx.find(it->idx)->second));
           fstpst.write(code_file, instruction_number);
         }
       // Insert current node into tt2
-      tt2.insert(*it);
+      tt2.insert(it);
     }
 }
 
@@ -1480,21 +1463,20 @@ ModelTree::writeJsonModelLocalVariables(ostream &output, deriv_node_temp_terms_t
   // Use an empty set for the temporary terms
   const temporary_terms_t tt;
 
-  for (size_t i = 0; i < equations.size(); i++)
-    equations[i]->collectVariables(eModelLocalVariable, used_local_vars);
+  for (auto equation : equations)
+    equation->collectVariables(eModelLocalVariable, used_local_vars);
 
   output << "\"model_local_variables\": [";
   bool printed = false;
-  for (vector<int>::const_iterator it = local_variables_vector.begin();
-       it != local_variables_vector.end(); it++)
-    if (used_local_vars.find(*it) != used_local_vars.end())
+  for (int it : local_variables_vector)
+    if (used_local_vars.find(it) != used_local_vars.end())
       {
         if (printed)
           output << ", ";
         else
           printed = true;
 
-        int id = *it;
+        int id = it;
         vector<string> efout;
         expr_t value = local_variables_table.find(id)->second;
         value->writeJsonExternalFunctionOutput(efout, tt, tef_terms);
@@ -1642,12 +1624,12 @@ ModelTree::Write_Inf_To_Bin_File(const string &basename,
       exit(EXIT_FAILURE);
     }
   u_count_int = 0;
-  for (first_derivatives_t::const_iterator it = first_derivatives.begin(); it != first_derivatives.end(); it++)
+  for (const auto & first_derivative : first_derivatives)
     {
-      int deriv_id = it->first.second;
+      int deriv_id = first_derivative.first.second;
       if (getTypeByDerivID(deriv_id) == eEndogenous)
         {
-          int eq = it->first.first;
+          int eq = first_derivative.first.first;
           int symb = getSymbIDByDerivID(deriv_id);
           int var = symbol_table.getTypeSpecificID(symb);
           int lag = getLagByDerivID(deriv_id);
@@ -1699,10 +1681,8 @@ ModelTree::writeLatexModelFile(const string &basename, ExprNodeOutputType output
          << "\\footnotesize" << endl;
 
   // Write model local variables
-  for (vector<int>::const_iterator it = local_variables_vector.begin();
-       it != local_variables_vector.end(); it++)
+  for (int id : local_variables_vector)
     {
-      int id = *it;
       expr_t value = local_variables_table.find(id)->second;
 
       content_output << "\\begin{dmath*}" << endl
@@ -1718,19 +1698,18 @@ ModelTree::writeLatexModelFile(const string &basename, ExprNodeOutputType output
       bool wrote_eq_tag = false;
       if (write_equation_tags)
         {
-          for (vector<pair<int, pair<string, string> > >::const_iterator iteqt = equation_tags.begin();
-               iteqt != equation_tags.end(); iteqt++)
-            if (iteqt->first == eq)
+          for (const auto & equation_tag : equation_tags)
+            if (equation_tag.first == eq)
               {
                 if (!wrote_eq_tag)
                   content_output << "\\noindent[";
                 else
                   content_output << ", ";
 
-                content_output << iteqt->second.first;
+                content_output << equation_tag.second.first;
 
-                if (!(iteqt->second.second.empty()))
-                  content_output << "= `" << iteqt->second.second << "'";
+                if (!(equation_tag.second.second.empty()))
+                  content_output << "= `" << equation_tag.second.second << "'";
 
                 wrote_eq_tag = true;
               }
@@ -1765,8 +1744,8 @@ void
 ModelTree::addEquation(expr_t eq, int lineno, const vector<pair<string, string> > &eq_tags)
 {
   int n = equations.size();
-  for (size_t i = 0; i < eq_tags.size(); i++)
-    equation_tags.push_back(make_pair(n, eq_tags[i]));
+  for (const auto & eq_tag : eq_tags)
+    equation_tags.push_back(make_pair(n, eq_tag));
   addEquation(eq, lineno);
 }
 
@@ -1853,11 +1832,8 @@ ModelTree::computeParamsDerivatives(int paramsDerivsOrder)
   set<int> deriv_id_set;
   addAllParamDerivId(deriv_id_set);
 
-  for (set<int>::const_iterator it = deriv_id_set.begin();
-       it != deriv_id_set.end(); it++)
+  for (int param : deriv_id_set)
     {
-      const int param = *it;
-
       for (int eq = 0; eq < (int) equations.size(); eq++)
         {
           expr_t d1 = equations[eq]->getDerivative(param);
@@ -1938,15 +1914,13 @@ ModelTree::computeParamsDerivativesTemporaryTerms()
   temp_terms_map[eJacobianParamsSecondDeriv] = params_derivs_temporary_terms_g12;
   temp_terms_map[eHessianParamsDeriv] = params_derivs_temporary_terms_g2;
 
-  for (first_derivatives_t::iterator it = residuals_params_derivatives.begin();
-       it != residuals_params_derivatives.end(); it++)
-    it->second->computeTemporaryTerms(reference_count,
+  for (auto & residuals_params_derivative : residuals_params_derivatives)
+    residuals_params_derivative.second->computeTemporaryTerms(reference_count,
                                       temp_terms_map,
                                       true, eResidualsParamsDeriv);
 
-  for (second_derivatives_t::iterator it = jacobian_params_derivatives.begin();
-       it != jacobian_params_derivatives.end(); it++)
-    it->second->computeTemporaryTerms(reference_count,
+  for (auto & jacobian_params_derivative : jacobian_params_derivatives)
+    jacobian_params_derivative.second->computeTemporaryTerms(reference_count,
                                       temp_terms_map,
                                       true, eJacobianParamsDeriv);
 
@@ -2054,10 +2028,9 @@ ModelTree::writeJsonModelEquations(ostream &output, bool residuals) const
           output << "\""
                  << ", \"line\": " << equations_lineno[eq];
 
-          for (vector<pair<int, pair<string, string> > >::const_iterator it = equation_tags.begin();
-               it != equation_tags.end(); it++)
-            if (it->first == eq)
-              eqtags.push_back(it->second);
+          for (const auto & equation_tag : equation_tags)
+            if (equation_tag.first == eq)
+              eqtags.push_back(equation_tag.second);
 
           if (!eqtags.empty())
             {
diff --git a/src/NumericalInitialization.cc b/src/NumericalInitialization.cc
index 4601634b..b3b39c98 100644
--- a/src/NumericalInitialization.cc
+++ b/src/NumericalInitialization.cc
@@ -106,12 +106,11 @@ InitOrEndValStatement::InitOrEndValStatement(const init_values_t &init_values_ar
 void
 InitOrEndValStatement::fillEvalContext(eval_context_t &eval_context) const
 {
-  for (init_values_t::const_iterator it = init_values.begin();
-       it != init_values.end(); it++)
+  for (const auto & init_value : init_values)
     {
       try
         {
-          eval_context[it->first] = (it->second)->eval(eval_context);
+          eval_context[init_value.first] = (init_value.second)->eval(eval_context);
         }
       catch (ExprNode::EvalException &e)
         {
@@ -138,10 +137,9 @@ InitOrEndValStatement::getUninitializedVariables(SymbolType type)
     }
 
   set<int>::iterator sit;
-  for (init_values_t::const_iterator it = init_values.begin();
-       it != init_values.end(); it++)
+  for (const auto & init_value : init_values)
     {
-      sit = unused.find(it->first);
+      sit = unused.find(init_value.first);
       if (sit != unused.end())
         unused.erase(sit);
     }
@@ -151,11 +149,10 @@ InitOrEndValStatement::getUninitializedVariables(SymbolType type)
 void
 InitOrEndValStatement::writeInitValues(ostream &output) const
 {
-  for (init_values_t::const_iterator it = init_values.begin();
-       it != init_values.end(); it++)
+  for (const auto & init_value : init_values)
     {
-      const int symb_id = it->first;
-      const expr_t expression = it->second;
+      const int symb_id = init_value.first;
+      const expr_t expression = init_value.second;
 
       SymbolType type = symbol_table.getType(symb_id);
       int tsid = symbol_table.getTypeSpecificID(symb_id) + 1;
@@ -203,16 +200,16 @@ InitValStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolidat
   if (endogs.size() > 0)
     {
       cerr << "ERROR: You have not set the following endogenous variables in initval:";
-      for (set<int>::const_iterator it = endogs.begin(); it != endogs.end(); it++)
-        cerr << " " << symbol_table.getName(*it);
+      for (int endog : endogs)
+        cerr << " " << symbol_table.getName(endog);
       cerr << endl;
     }
 
   if (exogs.size() > 0)
     {
       cerr << "ERROR: You have not set the following exogenous variables in initval:";
-      for (set<int>::const_iterator it = exogs.begin(); it != exogs.end(); it++)
-        cerr << " " << symbol_table.getName(*it);
+      for (int exog : exogs)
+        cerr << " " << symbol_table.getName(exog);
       cerr << endl;
     }
 
@@ -267,16 +264,16 @@ EndValStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolidati
   if (endogs.size() > 0)
     {
       cerr << "ERROR: You have not set the following endogenous variables in endval:";
-      for (set<int>::const_iterator it = endogs.begin(); it != endogs.end(); it++)
-        cerr << " " << symbol_table.getName(*it);
+      for (int endog : endogs)
+        cerr << " " << symbol_table.getName(endog);
       cerr << endl;
     }
 
   if (exogs.size() > 0)
     {
       cerr << "ERROR: You have not set the following exogenous variables in endval:";
-      for (set<int>::const_iterator it = exogs.begin(); it != exogs.end(); it++)
-        cerr << " " << symbol_table.getName(*it);
+      for (int exog : exogs)
+        cerr << " " << symbol_table.getName(exog);
       cerr << endl;
     }
 
@@ -325,14 +322,13 @@ HistValStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolidat
       set<int> unused_exo = symbol_table.getExogenous();
 
       set<int>::iterator sit;
-      for (hist_values_t::const_iterator it = hist_values.begin();
-           it != hist_values.end(); it++)
+      for (const auto & hist_value : hist_values)
         {
-          sit = unused_endo.find(it->first.first);
+          sit = unused_endo.find(hist_value.first.first);
           if (sit != unused_endo.end())
             unused_endo.erase(sit);
 
-          sit = unused_exo.find(it->first.first);
+          sit = unused_exo.find(hist_value.first.first);
           if (sit != unused_exo.end())
             unused_exo.erase(sit);
         }
@@ -340,16 +336,16 @@ HistValStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolidat
       if (unused_endo.size() > 0)
         {
           cerr << "ERROR: You have not set the following endogenous variables in histval:";
-          for (set<int>::const_iterator it = unused_endo.begin(); it != unused_endo.end(); it++)
-            cerr << " " << symbol_table.getName(*it);
+          for (int it : unused_endo)
+            cerr << " " << symbol_table.getName(it);
           cerr << endl;
         }
 
       if (unused_exo.size() > 0)
         {
           cerr << "ERROR: You have not set the following exogenous variables in endval:";
-          for (set<int>::const_iterator it = unused_exo.begin(); it != unused_exo.end(); it++)
-            cerr << " " << symbol_table.getName(*it);
+          for (int it : unused_exo)
+            cerr << " " << symbol_table.getName(it);
           cerr << endl;
         }
 
@@ -369,12 +365,11 @@ HistValStatement::writeOutput(ostream &output, const string &basename, bool mini
          << "M_.exo_histval = zeros(M_.exo_nbr,M_.maximum_lag);" << endl
          << "M_.exo_det_histval = zeros(M_.exo_det_nbr,M_.maximum_lag);" << endl;
 
-  for (hist_values_t::const_iterator it = hist_values.begin();
-       it != hist_values.end(); it++)
+  for (const auto & hist_value : hist_values)
     {
-      int symb_id = it->first.first;
-      int lag = it->first.second;
-      const expr_t expression = it->second;
+      int symb_id = hist_value.first.first;
+      int lag = hist_value.first.second;
+      const expr_t expression = hist_value.second;
 
       SymbolType type = symbol_table.getType(symb_id);
 
@@ -489,12 +484,11 @@ HomotopyStatement::writeOutput(ostream &output, const string &basename, bool min
          << "%" << endl
          << "options_.homotopy_values = [];" << endl;
 
-  for (homotopy_values_t::const_iterator it = homotopy_values.begin();
-       it != homotopy_values.end(); it++)
+  for (const auto & homotopy_value : homotopy_values)
     {
-      const int &symb_id = it->first;
-      const expr_t expression1 = it->second.first;
-      const expr_t expression2 = it->second.second;
+      const int &symb_id = homotopy_value.first;
+      const expr_t expression1 = homotopy_value.second.first;
+      const expr_t expression2 = homotopy_value.second.second;
 
       const SymbolType type = symbol_table.getType(symb_id);
       const int tsid = symbol_table.getTypeSpecificID(symb_id) + 1;
@@ -591,10 +585,9 @@ LoadParamsAndSteadyStateStatement::LoadParamsAndSteadyStateStatement(const strin
 void
 LoadParamsAndSteadyStateStatement::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const
 {
-  for (map<int, string>::const_iterator it = content.begin();
-       it != content.end(); it++)
+  for (const auto & it : content)
     {
-      switch (symbol_table.getType(it->first))
+      switch (symbol_table.getType(it.first))
         {
         case eParameter:
           output << "M_.params";
@@ -609,12 +602,12 @@ LoadParamsAndSteadyStateStatement::writeOutput(ostream &output, const string &ba
           output << "oo_.exo_det_steady_state";
           break;
         default:
-          cerr << "ERROR: Unsupported variable type for " << symbol_table.getName(it->first) << " in load_params_and_steady_state" << endl;
+          cerr << "ERROR: Unsupported variable type for " << symbol_table.getName(it.first) << " in load_params_and_steady_state" << endl;
           exit(EXIT_FAILURE);
         }
 
-      int tsid = symbol_table.getTypeSpecificID(it->first) + 1;
-      output << "(" << tsid << ") = " << it->second << ";" << endl;
+      int tsid = symbol_table.getTypeSpecificID(it.first) + 1;
+      output << "(" << tsid << ") = " << it.second << ";" << endl;
     }
 }
 
@@ -638,7 +631,6 @@ LoadParamsAndSteadyStateStatement::writeJsonOutput(ostream &output) const
 void
 LoadParamsAndSteadyStateStatement::fillEvalContext(eval_context_t &eval_context) const
 {
-  for (map<int, string>::const_iterator it = content.begin();
-       it != content.end(); it++)
-    eval_context[it->first] = atof(it->second.c_str());
+  for (const auto & it : content)
+    eval_context[it.first] = atof(it.second.c_str());
 }
diff --git a/src/ParsingDriver.cc b/src/ParsingDriver.cc
index a6e6bf3b..83a52c8b 100644
--- a/src/ParsingDriver.cc
+++ b/src/ParsingDriver.cc
@@ -215,12 +215,11 @@ ParsingDriver::declare_endogenous(string *name, string *tex_name, vector<pair<st
     delete tex_name;
   if (partition_value != NULL)
     {
-      for (vector<pair<string *, string *> *>::iterator it = partition_value->begin();
-           it != partition_value->end(); ++it)
+      for (auto & it : *partition_value)
         {
-          delete (*it)->first;
-          delete (*it)->second;
-          delete (*it);
+          delete it->first;
+          delete it->second;
+          delete it;
         }
       delete partition_value;
     }
@@ -252,12 +251,11 @@ ParsingDriver::declare_exogenous(string *name, string *tex_name, vector<pair<str
     delete tex_name;
   if (partition_value != NULL)
     {
-      for (vector<pair<string *, string *> *>::iterator it = partition_value->begin();
-           it != partition_value->end(); ++it)
+      for (auto & it : *partition_value)
         {
-          delete (*it)->first;
-          delete (*it)->second;
-          delete (*it);
+          delete it->first;
+          delete it->second;
+          delete it;
         }
       delete partition_value;
     }
@@ -272,12 +270,11 @@ ParsingDriver::declare_exogenous_det(string *name, string *tex_name, vector<pair
     delete tex_name;
   if (partition_value != NULL)
     {
-      for (vector<pair<string *, string *> *>::iterator it = partition_value->begin();
-           it != partition_value->end(); ++it)
+      for (auto & it : *partition_value)
         {
-          delete (*it)->first;
-          delete (*it)->second;
-          delete (*it);
+          delete it->first;
+          delete it->second;
+          delete it;
         }
       delete partition_value;
     }
@@ -292,12 +289,11 @@ ParsingDriver::declare_parameter(string *name, string *tex_name, vector<pair<str
     delete tex_name;
   if (partition_value != NULL)
     {
-      for (vector<pair<string *, string *> *>::iterator it = partition_value->begin();
-           it != partition_value->end(); ++it)
+      for (auto & it : *partition_value)
         {
-          delete (*it)->first;
-          delete (*it)->second;
-          delete (*it);
+          delete it->first;
+          delete it->second;
+          delete it;
         }
       delete partition_value;
     }
@@ -552,8 +548,8 @@ ParsingDriver::end_nonstationary_var(bool log_deflator, expr_t deflator)
 
   set<int> r;
   deflator->collectVariables(eEndogenous, r);
-  for (set<int>::const_iterator it = r.begin(); it != r.end(); ++it)
-    if (dynamic_model->isNonstationary(*it))
+  for (int it : r)
+    if (dynamic_model->isNonstationary(it))
       error("The deflator contains a non-stationary endogenous variable. This is not allowed. Please use only stationary endogenous and/or {log_}trend_vars.");
 
   declared_nonstationary_vars.clear();
@@ -2306,9 +2302,8 @@ ParsingDriver::run_identification()
 void
 ParsingDriver::add_mc_filename(string *filename, string *prior)
 {
-  for (ModelComparisonStatement::filename_list_t::iterator it = filename_list.begin();
-       it != filename_list.end(); it++)
-    if ((*it).first == *filename)
+  for (auto & it : filename_list)
+    if (it.first == *filename)
       error("model_comparison: filename " + *filename + " declared twice");
   filename_list.push_back(make_pair(*filename, *prior));
   delete filename;
@@ -2505,8 +2500,8 @@ ParsingDriver::svar()
 
   itv = options_list.vector_int_options.find("ms.equations");
   if (itv != options_list.vector_int_options.end())
-    for (vector<int>::const_iterator viit = itv->second.begin(); viit != itv->second.end(); viit++)
-      if (*viit <= 0)
+    for (int viit : itv->second)
+      if (viit <= 0)
         error("The value(s) passed to the equation option must be greater than zero.");
 
   mod_file->addStatement(new SvarStatement(options_list));
@@ -2688,27 +2683,26 @@ ParsingDriver::declare_and_init_model_local_variable(string *name, expr_t rhs)
 void
 ParsingDriver::change_type(SymbolType new_type, vector<string *> *var_list)
 {
-  for (vector<string *>::iterator it = var_list->begin();
-       it != var_list->end(); it++)
+  for (auto & it : *var_list)
     {
       int id;
       try
         {
-          id = mod_file->symbol_table.getID(**it);
+          id = mod_file->symbol_table.getID(*it);
         }
       catch (SymbolTable::UnknownSymbolNameException &e)
         {
-          error("Unknown variable " + **it);
+          error("Unknown variable " + *it);
         }
 
       // Check if symbol already used in a VariableNode
       if (mod_file->expressions_tree.isSymbolUsed(id)
           || mod_file->dynamic_model.isSymbolUsed(id))
-        error("You cannot modify the type of symbol " + **it + " after having used it in an expression");
+        error("You cannot modify the type of symbol " + *it + " after having used it in an expression");
 
       mod_file->symbol_table.changeType(id, new_type);
 
-      delete *it;
+      delete it;
     }
   delete var_list;
 }
@@ -3318,21 +3312,21 @@ ParsingDriver::add_steady_state_model_equal_multiple(expr_t expr)
   const vector<string> &symbs = symbol_list.get_symbols();
   vector<int> ids;
 
-  for (size_t i = 0; i < symbs.size(); i++)
+  for (const auto & symb : symbs)
     {
       int id;
       try
         {
-          id = mod_file->symbol_table.getID(symbs[i]);
+          id = mod_file->symbol_table.getID(symb);
         }
       catch (SymbolTable::UnknownSymbolNameException &e)
         {
           // Unknown symbol, declare it as a ModFileLocalVariable
-          id = mod_file->symbol_table.addSymbol(symbs[i], eModFileLocalVariable);
+          id = mod_file->symbol_table.addSymbol(symb, eModFileLocalVariable);
         }
       SymbolType type = mod_file->symbol_table.getType(id);
       if (type != eEndogenous && type != eModFileLocalVariable && type != eParameter)
-        error(symbs[i] + " has incorrect type");
+        error(symb + " has incorrect type");
       ids.push_back(id);
     }
 
diff --git a/src/Shocks.cc b/src/Shocks.cc
index cc8343c0..6f144bf6 100644
--- a/src/Shocks.cc
+++ b/src/Shocks.cc
@@ -39,17 +39,16 @@ AbstractShocksStatement::writeDetShocks(ostream &output) const
 {
   int exo_det_length = 0;
 
-  for (det_shocks_t::const_iterator it = det_shocks.begin();
-       it != det_shocks.end(); it++)
+  for (const auto & det_shock : det_shocks)
     {
-      int id = symbol_table.getTypeSpecificID(it->first) + 1;
-      bool exo_det = (symbol_table.getType(it->first) == eExogenousDet);
+      int id = symbol_table.getTypeSpecificID(det_shock.first) + 1;
+      bool exo_det = (symbol_table.getType(det_shock.first) == eExogenousDet);
 
-      for (size_t i = 0; i < it->second.size(); i++)
+      for (size_t i = 0; i < det_shock.second.size(); i++)
         {
-          const int &period1 = it->second[i].period1;
-          const int &period2 = it->second[i].period2;
-          const expr_t value = it->second[i].value;
+          const int &period1 = det_shock.second[i].period1;
+          const int &period2 = det_shock.second[i].period2;
+          const expr_t value = det_shock.second[i].value;
 
           output << "M_.det_shocks = [ M_.det_shocks;" << endl
                  << "struct('exo_det'," << (int) exo_det
@@ -317,35 +316,32 @@ ShocksStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolidati
   /* Error out if variables are not of the right type. This must be done here
      and not at parsing time (see #448).
      Also Determine if there is a calibrated measurement error */
-  for (var_and_std_shocks_t::const_iterator it = var_shocks.begin();
-       it != var_shocks.end(); it++)
+  for (auto var_shock : var_shocks)
     {
-      if (symbol_table.getType(it->first) != eExogenous
-          && !symbol_table.isObservedVariable(it->first))
+      if (symbol_table.getType(var_shock.first) != eExogenous
+          && !symbol_table.isObservedVariable(var_shock.first))
         {
           cerr << "shocks: setting a variance on '"
-               << symbol_table.getName(it->first) << "' is not allowed, because it is neither an exogenous variable nor an observed endogenous variable" << endl;
+               << symbol_table.getName(var_shock.first) << "' is not allowed, because it is neither an exogenous variable nor an observed endogenous variable" << endl;
           exit(EXIT_FAILURE);
         }
     }
 
-  for (var_and_std_shocks_t::const_iterator it = std_shocks.begin();
-       it != std_shocks.end(); it++)
+  for (auto std_shock : std_shocks)
     {
-      if (symbol_table.getType(it->first) != eExogenous
-          && !symbol_table.isObservedVariable(it->first))
+      if (symbol_table.getType(std_shock.first) != eExogenous
+          && !symbol_table.isObservedVariable(std_shock.first))
         {
           cerr << "shocks: setting a standard error on '"
-               << symbol_table.getName(it->first) << "' is not allowed, because it is neither an exogenous variable nor an observed endogenous variable" << endl;
+               << symbol_table.getName(std_shock.first) << "' is not allowed, because it is neither an exogenous variable nor an observed endogenous variable" << endl;
           exit(EXIT_FAILURE);
         }
     }
 
-  for (covar_and_corr_shocks_t::const_iterator it = covar_shocks.begin();
-       it != covar_shocks.end(); it++)
+  for (const auto & covar_shock : covar_shocks)
     {
-      int symb_id1 = it->first.first;
-      int symb_id2 = it->first.second;
+      int symb_id1 = covar_shock.first.first;
+      int symb_id2 = covar_shock.first.second;
 
       if (!((symbol_table.getType(symb_id1) == eExogenous
              && symbol_table.getType(symb_id2) == eExogenous)
@@ -359,11 +355,10 @@ ShocksStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolidati
         }
     }
 
-  for (covar_and_corr_shocks_t::const_iterator it = corr_shocks.begin();
-       it != corr_shocks.end(); it++)
+  for (const auto & corr_shock : corr_shocks)
     {
-      int symb_id1 = it->first.first;
-      int symb_id2 = it->first.second;
+      int symb_id1 = corr_shock.first.first;
+      int symb_id2 = corr_shock.first.second;
 
       if (!((symbol_table.getType(symb_id1) == eExogenous
              && symbol_table.getType(symb_id2) == eExogenous)
@@ -381,44 +376,36 @@ ShocksStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolidati
   mod_file_struct.calibrated_measurement_errors |= has_calibrated_measurement_errors();
 
   // Fill in mod_file_struct.parameters_with_shocks_values (related to #469)
-  for (var_and_std_shocks_t::const_iterator it = var_shocks.begin();
-       it != var_shocks.end(); ++it)
-    it->second->collectVariables(eParameter, mod_file_struct.parameters_within_shocks_values);
-  for (var_and_std_shocks_t::const_iterator it = std_shocks.begin();
-       it != std_shocks.end(); ++it)
-    it->second->collectVariables(eParameter, mod_file_struct.parameters_within_shocks_values);
-  for (covar_and_corr_shocks_t::const_iterator it = covar_shocks.begin();
-       it != covar_shocks.end(); ++it)
-    it->second->collectVariables(eParameter, mod_file_struct.parameters_within_shocks_values);
-  for (covar_and_corr_shocks_t::const_iterator it = corr_shocks.begin();
-       it != corr_shocks.end(); ++it)
-    it->second->collectVariables(eParameter, mod_file_struct.parameters_within_shocks_values);
+  for (auto var_shock : var_shocks)
+    var_shock.second->collectVariables(eParameter, mod_file_struct.parameters_within_shocks_values);
+  for (auto std_shock : std_shocks)
+    std_shock.second->collectVariables(eParameter, mod_file_struct.parameters_within_shocks_values);
+  for (const auto & covar_shock : covar_shocks)
+    covar_shock.second->collectVariables(eParameter, mod_file_struct.parameters_within_shocks_values);
+  for (const auto & corr_shock : corr_shocks)
+    corr_shock.second->collectVariables(eParameter, mod_file_struct.parameters_within_shocks_values);
 
 }
 
 bool
 ShocksStatement::has_calibrated_measurement_errors() const
 {
-  for (var_and_std_shocks_t::const_iterator it = var_shocks.begin();
-       it != var_shocks.end(); it++)
-    if (symbol_table.isObservedVariable(it->first))
+  for (auto var_shock : var_shocks)
+    if (symbol_table.isObservedVariable(var_shock.first))
       return true;
 
-  for (var_and_std_shocks_t::const_iterator it = std_shocks.begin();
-       it != std_shocks.end(); it++)
-    if (symbol_table.isObservedVariable(it->first))
+  for (auto std_shock : std_shocks)
+    if (symbol_table.isObservedVariable(std_shock.first))
       return true;
 
-  for (covar_and_corr_shocks_t::const_iterator it = covar_shocks.begin();
-       it != covar_shocks.end(); it++)
-    if (symbol_table.isObservedVariable(it->first.first)
-        || symbol_table.isObservedVariable(it->first.second))
+  for (const auto & covar_shock : covar_shocks)
+    if (symbol_table.isObservedVariable(covar_shock.first.first)
+        || symbol_table.isObservedVariable(covar_shock.first.second))
       return true;
 
-  for (covar_and_corr_shocks_t::const_iterator it = corr_shocks.begin();
-       it != corr_shocks.end(); it++)
-    if (symbol_table.isObservedVariable(it->first.first)
-        || symbol_table.isObservedVariable(it->first.second))
+  for (const auto & corr_shock : corr_shocks)
+    if (symbol_table.isObservedVariable(corr_shock.first.first)
+        || symbol_table.isObservedVariable(corr_shock.first.second))
       return true;
 
   return false;
@@ -455,14 +442,13 @@ ConditionalForecastPathsStatement::ConditionalForecastPathsStatement(const Abstr
 void
 ConditionalForecastPathsStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings)
 {
-  for (AbstractShocksStatement::det_shocks_t::const_iterator it = paths.begin();
-       it != paths.end(); it++)
+  for (const auto & path : paths)
     {
       int this_path_length = 0;
-      const vector<AbstractShocksStatement::DetShockElement> &elems = it->second;
-      for (int i = 0; i < (int) elems.size(); i++)
+      const vector<AbstractShocksStatement::DetShockElement> &elems = path.second;
+      for (auto elem : elems)
         // Period1 < Period2, as enforced in ParsingDriver::add_period()
-        this_path_length = max(this_path_length, elems[i].period2);
+        this_path_length = max(this_path_length, elem.period2);
       if (path_length == -1)
         path_length = this_path_length;
       else if (path_length != this_path_length)
@@ -490,11 +476,11 @@ ConditionalForecastPathsStatement::writeOutput(ostream &output, const string &ba
         output << "constrained_vars_ = [constrained_vars_; " << symbol_table.getTypeSpecificID(it->first) + 1 << "];" << endl;
 
       const vector<AbstractShocksStatement::DetShockElement> &elems = it->second;
-      for (int i = 0; i < (int) elems.size(); i++)
-        for (int j = elems[i].period1; j <= elems[i].period2; j++)
+      for (auto elem : elems)
+        for (int j = elem.period1; j <= elem.period2; j++)
           {
             output << "constrained_paths_(" << k << "," << j << ")=";
-            elems[i].value->writeOutput(output);
+            elem.value->writeOutput(output);
             output << ";" << endl;
           }
     }
@@ -510,9 +496,8 @@ void
 MomentCalibration::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const
 {
   output << "options_.endogenous_prior_restrictions.moment = {" << endl;
-  for (size_t i = 0; i < constraints.size(); i++)
+  for (const auto & c : constraints)
     {
-      const Constraint &c = constraints[i];
       output << "'" << symbol_table.getName(c.endo1) << "', "
              << "'" << symbol_table.getName(c.endo2) << "', "
              << c.lags << ", "
@@ -555,9 +540,8 @@ IrfCalibration::writeOutput(ostream &output, const string &basename, bool minima
   options_list.writeOutput(output);
 
   output << "options_.endogenous_prior_restrictions.irf = {" << endl;
-  for (size_t i = 0; i < constraints.size(); i++)
+  for (const auto & c : constraints)
     {
-      const Constraint &c = constraints[i];
       output << "'" << symbol_table.getName(c.endo) << "', "
              << "'" << symbol_table.getName(c.exo) << "', "
              << c.periods << ", "
@@ -620,8 +604,8 @@ ShockGroupsStatement::writeOutput(ostream &output, const string &basename, bool
                  << ".group" << i << ".label = '" << it->name << "';" << endl
                  << "M_.shock_groups." << name
                  << ".group" << i << ".shocks = {";
-          for (vector<string>::const_iterator it1 = it->list.begin(); it1 != it->list.end(); it1++)
-            output << " '" << *it1 << "'";
+          for (const auto & it1 : it->list)
+            output << " '" << it1 << "'";
           output << "};" << endl;
           i++;
         }
diff --git a/src/Statement.cc b/src/Statement.cc
index 5884c685..7ad2b905 100644
--- a/src/Statement.cc
+++ b/src/Statement.cc
@@ -142,57 +142,50 @@ VerbatimStatement::writeJsonOutput(ostream &output) const
 void
 OptionsList::writeOutput(ostream &output) const
 {
-  for (num_options_t::const_iterator it = num_options.begin();
-       it != num_options.end(); it++)
-    output << "options_." << it->first << " = " << it->second << ";" << endl;
+  for (const auto & num_option : num_options)
+    output << "options_." << num_option.first << " = " << num_option.second << ";" << endl;
 
-  for (paired_num_options_t::const_iterator it = paired_num_options.begin();
-       it != paired_num_options.end(); it++)
-    output << "options_." << it->first << " = [" << it->second.first << "; "
-           << it->second.second << "];" << endl;
+  for (const auto & paired_num_option : paired_num_options)
+    output << "options_." << paired_num_option.first << " = [" << paired_num_option.second.first << "; "
+           << paired_num_option.second.second << "];" << endl;
 
-  for (string_options_t::const_iterator it = string_options.begin();
-       it != string_options.end(); it++)
-    output << "options_." << it->first << " = '" << it->second << "';" << endl;
+  for (const auto & string_option : string_options)
+    output << "options_." << string_option.first << " = '" << string_option.second << "';" << endl;
 
-  for (date_options_t::const_iterator it = date_options.begin();
-       it != date_options.end(); it++)
-    output << "options_." << it->first << " = " << it->second << ";" << endl;
+  for (const auto & date_option : date_options)
+    output << "options_." << date_option.first << " = " << date_option.second << ";" << endl;
 
-  for (symbol_list_options_t::const_iterator it = symbol_list_options.begin();
-       it != symbol_list_options.end(); it++)
-    it->second.writeOutput("options_." + it->first, output);
+  for (const auto & symbol_list_option : symbol_list_options)
+    symbol_list_option.second.writeOutput("options_." + symbol_list_option.first, output);
 
-  for (vec_int_options_t::const_iterator it = vector_int_options.begin();
-       it != vector_int_options.end(); it++)
+  for (const auto & vector_int_option : vector_int_options)
     {
-      output << "options_." << it->first << " = ";
-      if (it->second.size() > 1)
+      output << "options_." << vector_int_option.first << " = ";
+      if (vector_int_option.second.size() > 1)
         {
           output << "[";
-          for (vector<int>::const_iterator viit = it->second.begin();
-               viit != it->second.end(); viit++)
+          for (vector<int>::const_iterator viit = vector_int_option.second.begin();
+               viit != vector_int_option.second.end(); viit++)
             output << *viit << ";";
           output << "];" << endl;
         }
       else
-        output << it->second.front() << ";" << endl;
+        output << vector_int_option.second.front() << ";" << endl;
     }
 
-  for (vec_str_options_t::const_iterator it = vector_str_options.begin();
-       it != vector_str_options.end(); it++)
+  for (const auto & vector_str_option : vector_str_options)
     {
-      output << "options_." << it->first << " = ";
-      if (it->second.size() > 1)
+      output << "options_." << vector_str_option.first << " = ";
+      if (vector_str_option.second.size() > 1)
         {
           output << "{";
-          for (vector<string>::const_iterator viit = it->second.begin();
-               viit != it->second.end(); viit++)
+          for (vector<string>::const_iterator viit = vector_str_option.second.begin();
+               viit != vector_str_option.second.end(); viit++)
             output << "'" << *viit << "';";
           output << "};" << endl;
         }
       else
-        output << it->second.front() << ";" << endl;
+        output << vector_str_option.second.front() << ";" << endl;
     }
 }
 
@@ -210,57 +203,50 @@ OptionsList::writeOutput(ostream &output, const string &option_group) const
   else
     output << option_group << " = struct();" << endl;
 
-  for (num_options_t::const_iterator it = num_options.begin();
-       it != num_options.end(); it++)
-    output << option_group << "." << it->first << " = " << it->second << ";" << endl;
+  for (const auto & num_option : num_options)
+    output << option_group << "." << num_option.first << " = " << num_option.second << ";" << endl;
 
-  for (paired_num_options_t::const_iterator it = paired_num_options.begin();
-       it != paired_num_options.end(); it++)
-    output << option_group << "." << it->first << " = [" << it->second.first << "; "
-           << it->second.second << "];" << endl;
+  for (const auto & paired_num_option : paired_num_options)
+    output << option_group << "." << paired_num_option.first << " = [" << paired_num_option.second.first << "; "
+           << paired_num_option.second.second << "];" << endl;
 
-  for (string_options_t::const_iterator it = string_options.begin();
-       it != string_options.end(); it++)
-    output << option_group << "." << it->first << " = '" << it->second << "';" << endl;
+  for (const auto & string_option : string_options)
+    output << option_group << "." << string_option.first << " = '" << string_option.second << "';" << endl;
 
-  for (date_options_t::const_iterator it = date_options.begin();
-       it != date_options.end(); it++)
-    output << option_group << "." << it->first << " = " << it->second << ";" << endl;
+  for (const auto & date_option : date_options)
+    output << option_group << "." << date_option.first << " = " << date_option.second << ";" << endl;
 
-  for (symbol_list_options_t::const_iterator it = symbol_list_options.begin();
-       it != symbol_list_options.end(); it++)
-    it->second.writeOutput(option_group + "." + it->first, output);
+  for (const auto & symbol_list_option : symbol_list_options)
+    symbol_list_option.second.writeOutput(option_group + "." + symbol_list_option.first, output);
 
-  for (vec_int_options_t::const_iterator it = vector_int_options.begin();
-       it != vector_int_options.end(); it++)
+  for (const auto & vector_int_option : vector_int_options)
     {
-      output << option_group << "." << it->first << " = ";
-      if (it->second.size() > 1)
+      output << option_group << "." << vector_int_option.first << " = ";
+      if (vector_int_option.second.size() > 1)
         {
           output << "[";
-          for (vector<int>::const_iterator viit = it->second.begin();
-               viit != it->second.end(); viit++)
+          for (vector<int>::const_iterator viit = vector_int_option.second.begin();
+               viit != vector_int_option.second.end(); viit++)
             output << *viit << ";";
           output << "];" << endl;
         }
       else
-        output <<  it->second.front() << ";" << endl;
+        output <<  vector_int_option.second.front() << ";" << endl;
     }
 
-  for (vec_str_options_t::const_iterator it = vector_str_options.begin();
-       it != vector_str_options.end(); it++)
+  for (const auto & vector_str_option : vector_str_options)
     {
-      output << option_group << "." << it->first << " = ";
-      if (it->second.size() > 1)
+      output << option_group << "." << vector_str_option.first << " = ";
+      if (vector_str_option.second.size() > 1)
         {
           output << "{";
-          for (vector<string>::const_iterator viit = it->second.begin();
-               viit != it->second.end(); viit++)
+          for (vector<string>::const_iterator viit = vector_str_option.second.begin();
+               viit != vector_str_option.second.end(); viit++)
             output << "'" << *viit << "';";
           output << "};" << endl;
         }
       else
-        output <<  it->second.front() << ";" << endl;
+        output <<  vector_str_option.second.front() << ";" << endl;
     }
 }
 
diff --git a/src/StaticModel.cc b/src/StaticModel.cc
index 46fb73e8..7186c4b0 100644
--- a/src/StaticModel.cc
+++ b/src/StaticModel.cc
@@ -188,9 +188,8 @@ StaticModel::computeTemporaryTermsMapping(temporary_terms_t &temporary_terms, ma
 {
   // Add a mapping form node ID to temporary terms order
   int j = 0;
-  for (temporary_terms_t::const_iterator it = temporary_terms.begin();
-       it != temporary_terms.end(); it++)
-    map_idx[(*it)->idx] = j++;
+  for (auto temporary_term : temporary_terms)
+    map_idx[temporary_term->idx] = j++;
 }
 
 void
@@ -264,9 +263,8 @@ StaticModel::writeModelEquationsOrdered_M(const string &static_basename) const
       if (v_temporary_terms_inuse[block].size())
         {
           tmp_output.str("");
-          for (temporary_terms_inuse_t::const_iterator it = v_temporary_terms_inuse[block].begin();
-               it != v_temporary_terms_inuse[block].end(); it++)
-            tmp_output << " T" << *it;
+          for (int it : v_temporary_terms_inuse[block])
+            tmp_output << " T" << it;
           output << "  global" << tmp_output.str() << ";\n";
         }
 
@@ -284,18 +282,17 @@ StaticModel::writeModelEquationsOrdered_M(const string &static_basename) const
           if (v_temporary_terms[block].size())
             {
               output << "  " << "% //Temporary variables" << endl;
-              for (temporary_terms_t::const_iterator it = v_temporary_terms[block][i].begin();
-                   it != v_temporary_terms[block][i].end(); it++)
+              for (auto it : v_temporary_terms[block][i])
                 {
-                  if (dynamic_cast<AbstractExternalFunctionNode *>(*it) != NULL)
-                    (*it)->writeExternalFunctionOutput(output, local_output_type, tt2, {}, tef_terms);
+                  if (dynamic_cast<AbstractExternalFunctionNode *>(it) != NULL)
+                    it->writeExternalFunctionOutput(output, local_output_type, tt2, {}, tef_terms);
 
                   output << "  " <<  sps;
-                  (*it)->writeOutput(output, local_output_type, local_temporary_terms, {}, tef_terms);
+                  it->writeOutput(output, local_output_type, local_temporary_terms, {}, tef_terms);
                   output << " = ";
-                  (*it)->writeOutput(output, local_output_type, tt2, {}, tef_terms);
+                  it->writeOutput(output, local_output_type, tt2, {}, tef_terms);
                   // Insert current node into tt2
-                  tt2.insert(*it);
+                  tt2.insert(it);
                   output << ";" << endl;
                 }
             }
@@ -442,9 +439,8 @@ StaticModel::writeModelEquationsCode(const string file_name, const string bin_ba
 
   // Add a mapping form node ID to temporary terms order
   int j = 0;
-  for (temporary_terms_t::const_iterator it = temporary_terms.begin();
-       it != temporary_terms.end(); it++)
-    map_idx[(*it)->idx] = j++;
+  for (auto temporary_term : temporary_terms)
+    map_idx[temporary_term->idx] = j++;
   compileTemporaryTerms(code_file, instruction_number, temporary_terms, map_idx, false, false);
 
   compileModelEquations(code_file, instruction_number, temporary_terms, map_idx, false, false);
@@ -461,14 +457,13 @@ StaticModel::writeModelEquationsCode(const string file_name, const string bin_ba
   vector<vector<pair<int, int> > > derivatives;
   derivatives.resize(symbol_table.endo_nbr());
   count_u = symbol_table.endo_nbr();
-  for (first_derivatives_t::const_iterator it = first_derivatives.begin();
-       it != first_derivatives.end(); it++)
+  for (const auto & first_derivative : first_derivatives)
     {
-      int deriv_id = it->first.second;
+      int deriv_id = first_derivative.first.second;
       if (getTypeByDerivID(deriv_id) == eEndogenous)
         {
-          expr_t d1 = it->second;
-          unsigned int eq = it->first.first;
+          expr_t d1 = first_derivative.second;
+          unsigned int eq = first_derivative.first.first;
           int symb = getSymbIDByDerivID(deriv_id);
           unsigned int var = symbol_table.getTypeSpecificID(symb);
           FNUMEXPR_ fnumexpr(FirstEndoDerivative, eq, var);
@@ -529,14 +524,13 @@ StaticModel::writeModelEquationsCode(const string file_name, const string bin_ba
   tt3.clear();
 
   // The Jacobian if we have to solve the block determinsitic bloc
-  for (first_derivatives_t::const_iterator it = first_derivatives.begin();
-       it != first_derivatives.end(); it++)
+  for (const auto & first_derivative : first_derivatives)
     {
-      int deriv_id = it->first.second;
+      int deriv_id = first_derivative.first.second;
       if (getTypeByDerivID(deriv_id) == eEndogenous)
         {
-          expr_t d1 = it->second;
-          unsigned int eq = it->first.first;
+          expr_t d1 = first_derivative.second;
+          unsigned int eq = first_derivative.first.first;
           int symb = getSymbIDByDerivID(deriv_id);
           unsigned int var = symbol_table.getTypeSpecificID(symb);
           FNUMEXPR_ fnumexpr(FirstEndoDerivative, eq, var);
@@ -656,19 +650,18 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
           tt2.clear();
           if (v_temporary_terms[block].size())
             {
-              for (temporary_terms_t::const_iterator it = v_temporary_terms[block][i].begin();
-                   it != v_temporary_terms[block][i].end(); it++)
+              for (auto it : v_temporary_terms[block][i])
                 {
-                  if (dynamic_cast<AbstractExternalFunctionNode *>(*it) != NULL)
-                    (*it)->compileExternalFunctionOutput(code_file, instruction_number, false, tt2, map_idx, false, false, tef_terms);
+                  if (dynamic_cast<AbstractExternalFunctionNode *>(it) != NULL)
+                    it->compileExternalFunctionOutput(code_file, instruction_number, false, tt2, map_idx, false, false, tef_terms);
 
-                  FNUMEXPR_ fnumexpr(TemporaryTerm, (int)(map_idx.find((*it)->idx)->second));
+                  FNUMEXPR_ fnumexpr(TemporaryTerm, (int)(map_idx.find(it->idx)->second));
                   fnumexpr.write(code_file, instruction_number);
-                  (*it)->compile(code_file, instruction_number, false, tt2, map_idx, false, false, tef_terms);
-                  FSTPST_ fstpst((int)(map_idx.find((*it)->idx)->second));
+                  it->compile(code_file, instruction_number, false, tt2, map_idx, false, false, tef_terms);
+                  FSTPST_ fstpst((int)(map_idx.find(it->idx)->second));
                   fstpst.write(code_file, instruction_number);
                   // Insert current node into tt2
-                  tt2.insert(*it);
+                  tt2.insert(it);
                 }
             }
 
@@ -1033,15 +1026,14 @@ map<pair<int, pair<int, int > >, expr_t>
 StaticModel::collect_first_order_derivatives_endogenous()
 {
   map<pair<int, pair<int, int > >, expr_t> endo_derivatives;
-  for (first_derivatives_t::iterator it2 = first_derivatives.begin();
-       it2 != first_derivatives.end(); it2++)
+  for (auto & first_derivative : first_derivatives)
     {
-      if (getTypeByDerivID(it2->first.second) == eEndogenous)
+      if (getTypeByDerivID(first_derivative.first.second) == eEndogenous)
         {
-          int eq = it2->first.first;
-          int var = symbol_table.getTypeSpecificID(getSymbIDByDerivID(it2->first.second));
+          int eq = first_derivative.first.first;
+          int var = symbol_table.getTypeSpecificID(getSymbIDByDerivID(first_derivative.first.second));
           int lag = 0;
-          endo_derivatives[make_pair(eq, make_pair(var, lag))] = it2->second;
+          endo_derivatives[make_pair(eq, make_pair(var, lag))] = first_derivative.second;
         }
     }
   return endo_derivatives;
@@ -1059,9 +1051,9 @@ StaticModel::computingPass(const eval_context_t &eval_context, bool no_tmp_terms
       neweqs.push_back(dynamic_cast<BinaryOpNode *>(eq_tmp->toStatic(*this)));
     }
 
-  for (unsigned int eq = 0; eq < aux_equations.size(); eq++)
+  for (auto & aux_equation : aux_equations)
     {
-      expr_t eq_tmp = aux_equations[eq]->substituteStaticAuxiliaryDefinition();
+      expr_t eq_tmp = aux_equation->substituteStaticAuxiliaryDefinition();
       neweqs.push_back(dynamic_cast<BinaryOpNode *>(eq_tmp->toStatic(*this)));
     }
 
@@ -1356,9 +1348,8 @@ StaticModel::writeStaticModel(const string &basename,
   deriv_node_temp_terms_t tef_terms;
   temporary_terms_t temp_term_union;
 
-  for (map<expr_t, expr_t>::const_iterator it = temporary_terms_mlv.begin();
-       it != temporary_terms_mlv.end(); it++)
-    temp_term_union.insert(it->first);
+  for (auto it : temporary_terms_mlv)
+    temp_term_union.insert(it.first);
   writeModelLocalVariableTemporaryTerms(temp_term_union, temporary_terms_mlv,
                                         model_tt_output, output_type, tef_terms);
 
@@ -1383,12 +1374,11 @@ StaticModel::writeStaticModel(const string &basename,
                           jacobian_tt_output, output_type, tef_terms);
       temp_term_union.insert(temporary_terms_g1.begin(), temporary_terms_g1.end());
     }
-  for (first_derivatives_t::const_iterator it = first_derivatives.begin();
-       it != first_derivatives.end(); it++)
+  for (const auto & first_derivative : first_derivatives)
     {
-      int eq = it->first.first;
-      int symb_id = getSymbIDByDerivID(it->first.second);
-      expr_t d1 = it->second;
+      int eq = first_derivative.first.first;
+      int symb_id = getSymbIDByDerivID(first_derivative.first.second);
+      expr_t d1 = first_derivative.second;
 
       jacobianHelper(jacobian_output, eq, symbol_table.getTypeSpecificID(symb_id), output_type);
       jacobian_output << "=";
@@ -1408,13 +1398,12 @@ StaticModel::writeStaticModel(const string &basename,
       temp_term_union.insert(temporary_terms_g2.begin(), temporary_terms_g2.end());
 
       int k = 0; // Keep the line of a 2nd derivative in v2
-      for (second_derivatives_t::const_iterator it = second_derivatives.begin();
-           it != second_derivatives.end(); it++)
+      for (const auto & second_derivative : second_derivatives)
         {
-          int eq = it->first.first;
-          int symb_id1 = getSymbIDByDerivID(it->first.second.first);
-          int symb_id2 = getSymbIDByDerivID(it->first.second.second);
-          expr_t d2 = it->second;
+          int eq = second_derivative.first.first;
+          int symb_id1 = getSymbIDByDerivID(second_derivative.first.second.first);
+          int symb_id2 = getSymbIDByDerivID(second_derivative.first.second.second);
+          expr_t d2 = second_derivative.second;
 
           int tsid1 = symbol_table.getTypeSpecificID(symb_id1);
           int tsid2 = symbol_table.getTypeSpecificID(symb_id2);
@@ -1478,14 +1467,13 @@ StaticModel::writeStaticModel(const string &basename,
       temp_term_union.insert(temporary_terms_g3.begin(), temporary_terms_g3.end());
 
       int k = 0; // Keep the line of a 3rd derivative in v3
-      for (third_derivatives_t::const_iterator it = third_derivatives.begin();
-           it != third_derivatives.end(); it++)
+      for (const auto & third_derivative : third_derivatives)
         {
-          int eq = it->first.first;
-          int var1 = it->first.second.first;
-          int var2 = it->first.second.second.first;
-          int var3 = it->first.second.second.second;
-          expr_t d3 = it->second;
+          int eq = third_derivative.first.first;
+          int var1 = third_derivative.first.second.first;
+          int var2 = third_derivative.first.second.second.first;
+          int var3 = third_derivative.first.second.second.second;
+          expr_t d3 = third_derivative.second;
 
           int id1 = getSymbIDByDerivID(var1);
           int id2 = getSymbIDByDerivID(var2);
@@ -1525,10 +1513,10 @@ StaticModel::writeStaticModel(const string &basename,
           cols.insert(id3 * hessianColsNbr + id2 * JacobianColsNbr + id1);
 
           int k2 = 1; // Keeps the offset of the permutation relative to k
-          for (set<int>::iterator it2 = cols.begin(); it2 != cols.end(); it2++)
-            if (*it2 != ref_col)
+          for (int col : cols)
+            if (col != ref_col)
               if (output_type == oJuliaStaticModel)
-                third_derivatives_output << "  @inbounds g3[" << eq + 1 << "," << *it2 + 1 << "] = "
+                third_derivatives_output << "  @inbounds g3[" << eq + 1 << "," << col + 1 << "] = "
                                          << for_sym.str() << endl;
               else
                 {
@@ -1536,7 +1524,7 @@ StaticModel::writeStaticModel(const string &basename,
                   third_derivatives_output << "=" << eq + 1 << ";" << endl;
 
                   sparseHelper(3, third_derivatives_output, k+k2, 1, output_type);
-                  third_derivatives_output << "=" << *it2 + 1 << ";" << endl;
+                  third_derivatives_output << "=" << col + 1 << ";" << endl;
 
                   sparseHelper(3, third_derivatives_output, k+k2, 2, output_type);
                   third_derivatives_output << "=";
@@ -2042,8 +2030,8 @@ StaticModel::writeStaticFile(const string &basename, bool block, bool bytecode,
 bool
 StaticModel::exoPresentInEqs() const
 {
-  for (int i = 0; i < (int) equations.size(); i++)
-    if (equations[i]->containsExogenous())
+  for (auto equation : equations)
+    if (equation->containsExogenous())
       return true;
   return false;
 }
@@ -2143,13 +2131,12 @@ StaticModel::writeOutput(ostream &output, bool block) const
   output << "];\n";
 
   map<pair<int, int>,  int>  row_incidence;
-  for (first_derivatives_t::const_iterator it = first_derivatives.begin();
-       it != first_derivatives.end(); it++)
+  for (const auto & first_derivative : first_derivatives)
     {
-      int deriv_id = it->first.second;
+      int deriv_id = first_derivative.first.second;
       if (getTypeByDerivID(deriv_id) == eEndogenous)
         {
-          int eq = it->first.first;
+          int eq = first_derivative.first.first;
           int symb = getSymbIDByDerivID(deriv_id);
           int var = symbol_table.getTypeSpecificID(symb);
           //int lag = getLagByDerivID(deriv_id);
@@ -2360,11 +2347,10 @@ StaticModel::collect_block_first_order_derivatives()
   derivative_endo = vector<derivative_t>(nb_blocks);
   endo_max_leadlag_block = vector<pair<int, int> >(nb_blocks, make_pair(0, 0));
   max_leadlag_block = vector<pair<int, int> >(nb_blocks, make_pair(0, 0));
-  for (first_derivatives_t::iterator it2 = first_derivatives.begin();
-       it2 != first_derivatives.end(); it2++)
+  for (auto & first_derivative : first_derivatives)
     {
-      int eq = it2->first.first;
-      int var = symbol_table.getTypeSpecificID(getSymbIDByDerivID(it2->first.second));
+      int eq = first_derivative.first.first;
+      int var = symbol_table.getTypeSpecificID(getSymbIDByDerivID(first_derivative.first.second));
       int lag = 0;
       int block_eq = equation_2_block[eq];
       int block_var = variable_2_block[var];
@@ -2374,7 +2360,7 @@ StaticModel::collect_block_first_order_derivatives()
       endo_max_leadlag_block[block_eq] = make_pair(0, 0);
       derivative_t tmp_derivative;
       lag_var_t lag_var;
-      if (getTypeByDerivID(it2->first.second) == eEndogenous && block_eq == block_var)
+      if (getTypeByDerivID(first_derivative.first.second) == eEndogenous && block_eq == block_var)
         {
           tmp_derivative = derivative_endo[block_eq];
           tmp_derivative[make_pair(lag, make_pair(eq, var))] = first_derivatives[make_pair(eq, getDerivID(symbol_table.getID(eEndogenous, var), lag))];
@@ -2392,9 +2378,9 @@ StaticModel::writeLatexFile(const string &basename, bool write_equation_tags) co
 void
 StaticModel::writeAuxVarInitval(ostream &output, ExprNodeOutputType output_type) const
 {
-  for (int i = 0; i < (int) aux_equations.size(); i++)
+  for (auto aux_equation : aux_equations)
     {
-      dynamic_cast<ExprNode *>(aux_equations[i])->writeOutput(output, output_type);
+      dynamic_cast<ExprNode *>(aux_equation)->writeOutput(output, output_type);
       output << ";" << endl;
     }
 }
@@ -2435,12 +2421,12 @@ void
 StaticModel::writeAuxVarRecursiveDefinitions(ostream &output, ExprNodeOutputType output_type) const
 {
   deriv_node_temp_terms_t tef_terms;
-  for (int i = 0; i < (int) aux_equations.size(); i++)
-    if (dynamic_cast<ExprNode *>(aux_equations[i])->containsExternalFunction())
-      dynamic_cast<ExprNode *>(aux_equations[i])->writeExternalFunctionOutput(output, oMatlabStaticModel, {}, {}, tef_terms);
-  for (int i = 0; i < (int) aux_equations.size(); i++)
+  for (auto aux_equation : aux_equations)
+    if (dynamic_cast<ExprNode *>(aux_equation)->containsExternalFunction())
+      dynamic_cast<ExprNode *>(aux_equation)->writeExternalFunctionOutput(output, oMatlabStaticModel, {}, {}, tef_terms);
+  for (auto aux_equation : aux_equations)
     {
-      dynamic_cast<ExprNode *>(aux_equations[i]->substituteStaticAuxiliaryDefinition())->writeOutput(output, output_type);
+      dynamic_cast<ExprNode *>(aux_equation->substituteStaticAuxiliaryDefinition())->writeOutput(output, output_type);
       output << ";" << endl;
     }
 }
@@ -2451,14 +2437,14 @@ StaticModel::writeLatexAuxVarRecursiveDefinitions(ostream &output) const
   deriv_node_temp_terms_t tef_terms;
   temporary_terms_t temporary_terms;
     temporary_terms_idxs_t temporary_terms_idxs;
-  for (int i = 0; i < (int) aux_equations.size(); i++)
-    if (dynamic_cast<ExprNode *>(aux_equations[i])->containsExternalFunction())
-      dynamic_cast<ExprNode *>(aux_equations[i])->writeExternalFunctionOutput(output, oLatexStaticModel,
+  for (auto aux_equation : aux_equations)
+    if (dynamic_cast<ExprNode *>(aux_equation)->containsExternalFunction())
+      dynamic_cast<ExprNode *>(aux_equation)->writeExternalFunctionOutput(output, oLatexStaticModel,
                                                                               temporary_terms, temporary_terms_idxs, tef_terms);
-  for (int i = 0; i < (int) aux_equations.size(); i++)
+  for (auto aux_equation : aux_equations)
     {
       output << "\\begin{dmath}" << endl;
-      dynamic_cast<ExprNode *>(aux_equations[i]->substituteStaticAuxiliaryDefinition())->writeOutput(output, oLatexStaticModel);
+      dynamic_cast<ExprNode *>(aux_equation->substituteStaticAuxiliaryDefinition())->writeOutput(output, oLatexStaticModel);
       output << endl << "\\end{dmath}" << endl;
     }
 }
@@ -2469,11 +2455,11 @@ StaticModel::writeJsonAuxVarRecursiveDefinitions(ostream &output) const
   deriv_node_temp_terms_t tef_terms;
   temporary_terms_t temporary_terms;
 
-  for (int i = 0; i < (int) aux_equations.size(); i++)
-    if (dynamic_cast<ExprNode *>(aux_equations[i])->containsExternalFunction())
+  for (auto aux_equation : aux_equations)
+    if (dynamic_cast<ExprNode *>(aux_equation)->containsExternalFunction())
       {
         vector<string> efout;
-        dynamic_cast<ExprNode *>(aux_equations[i])->writeJsonExternalFunctionOutput(efout,
+        dynamic_cast<ExprNode *>(aux_equation)->writeJsonExternalFunctionOutput(efout,
                                                                                     temporary_terms,
                                                                                     tef_terms,
                                                                                     false);
@@ -2485,12 +2471,12 @@ StaticModel::writeJsonAuxVarRecursiveDefinitions(ostream &output) const
           }
       }
 
-  for (int i = 0; i < (int) aux_equations.size(); i++)
+  for (auto aux_equation : aux_equations)
     {
       output << ", {\"lhs\": \"";
-      aux_equations[i]->get_arg1()->writeJsonOutput(output, temporary_terms, tef_terms, false);
+      aux_equation->get_arg1()->writeJsonOutput(output, temporary_terms, tef_terms, false);
       output << "\", \"rhs\": \"";
-      dynamic_cast<BinaryOpNode *>(aux_equations[i]->substituteStaticAuxiliaryDefinition())->get_arg2()->writeJsonOutput(output, temporary_terms, tef_terms, false);
+      dynamic_cast<BinaryOpNode *>(aux_equation->substituteStaticAuxiliaryDefinition())->get_arg2()->writeJsonOutput(output, temporary_terms, tef_terms, false);
       output << "\"}";
     }
 }
@@ -2519,12 +2505,11 @@ StaticModel::writeParamsDerivativesFile(const string &basename, bool julia) cons
 
   writeTemporaryTerms(params_derivs_temporary_terms, {}, params_derivs_temporary_terms_idxs, model_output, output_type, tef_terms);
 
-  for (first_derivatives_t::const_iterator it = residuals_params_derivatives.begin();
-       it != residuals_params_derivatives.end(); it++)
+  for (const auto & residuals_params_derivative : residuals_params_derivatives)
     {
-      int eq = it->first.first;
-      int param = it->first.second;
-      expr_t d1 = it->second;
+      int eq = residuals_params_derivative.first.first;
+      int param = residuals_params_derivative.first.second;
+      expr_t d1 = residuals_params_derivative.second;
 
       int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1;
 
@@ -2535,13 +2520,12 @@ StaticModel::writeParamsDerivativesFile(const string &basename, bool julia) cons
       jacobian_output << ";" << endl;
     }
 
-  for (second_derivatives_t::const_iterator it = jacobian_params_derivatives.begin();
-       it != jacobian_params_derivatives.end(); it++)
+  for (const auto & jacobian_params_derivative : jacobian_params_derivatives)
     {
-      int eq = it->first.first;
-      int var = it->first.second.first;
-      int param = it->first.second.second;
-      expr_t d2 = it->second;
+      int eq = jacobian_params_derivative.first.first;
+      int var = jacobian_params_derivative.first.second.first;
+      int param = jacobian_params_derivative.first.second.second;
+      expr_t d2 = jacobian_params_derivative.second;
 
       int var_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(var)) + 1;
       int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1;
diff --git a/src/SteadyStateModel.cc b/src/SteadyStateModel.cc
index 37a9dd34..24bb2d53 100644
--- a/src/SteadyStateModel.cc
+++ b/src/SteadyStateModel.cc
@@ -45,12 +45,12 @@ SteadyStateModel::addDefinition(int symb_id, expr_t expr)
 void
 SteadyStateModel::addMultipleDefinitions(const vector<int> &symb_ids, expr_t expr)
 {
-  for (size_t i = 0; i < symb_ids.size(); i++)
+  for (int symb_id : symb_ids)
     {
-      AddVariable(symb_ids[i]); // Create the variable nodes to be used in write method
-      assert(symbol_table.getType(symb_ids[i]) == eEndogenous
-             || symbol_table.getType(symb_ids[i]) == eModFileLocalVariable
-             || symbol_table.getType(symb_ids[i]) == eParameter);
+      AddVariable(symb_id); // Create the variable nodes to be used in write method
+      assert(symbol_table.getType(symb_id) == eEndogenous
+             || symbol_table.getType(symb_id) == eModFileLocalVariable
+             || symbol_table.getType(symb_id) == eParameter);
     }
   def_table.push_back(make_pair(symb_ids, expr));
 }
@@ -64,29 +64,28 @@ SteadyStateModel::checkPass(ModFileStructure &mod_file_struct, WarningConsolidat
   mod_file_struct.steady_state_model_present = true;
   vector<int> so_far_defined;
 
-  for (size_t i = 0; i < def_table.size(); i++)
+  for (const auto & i : def_table)
     {
-      const vector<int> &symb_ids = def_table[i].first;
+      const vector<int> &symb_ids = i.first;
 
       // Check that symbols are not already defined
-      for (size_t j = 0; j < symb_ids.size(); j++)
-        if (find(so_far_defined.begin(), so_far_defined.end(), symb_ids[j])
+      for (int symb_id : symb_ids)
+        if (find(so_far_defined.begin(), so_far_defined.end(), symb_id)
             != so_far_defined.end())
-          warnings << "WARNING: in the 'steady_state_model' block, variable '" << symbol_table.getName(symb_ids[j]) << "' is declared twice" << endl;
+          warnings << "WARNING: in the 'steady_state_model' block, variable '" << symbol_table.getName(symb_id) << "' is declared twice" << endl;
 
       // Check that expression has no undefined symbol
       if (!mod_file_struct.ramsey_model_present)
         {
           set<int> used_symbols;
-          const expr_t &expr = def_table[i].second;
+          const expr_t &expr = i.second;
           expr->collectVariables(eEndogenous, used_symbols);
           expr->collectVariables(eModFileLocalVariable, used_symbols);
-          for (set<int>::const_iterator it = used_symbols.begin();
-               it != used_symbols.end(); ++it)
-            if (find(so_far_defined.begin(), so_far_defined.end(), *it)
+          for (int used_symbol : used_symbols)
+            if (find(so_far_defined.begin(), so_far_defined.end(), used_symbol)
                 == so_far_defined.end())
               {
-                cerr << "ERROR: in the 'steady_state_model' block, variable '" << symbol_table.getName(*it)
+                cerr << "ERROR: in the 'steady_state_model' block, variable '" << symbol_table.getName(used_symbol)
                      << "' is undefined in the declaration of variable '" << symbol_table.getName(symb_ids[0]) << "'" << endl;
                 exit(EXIT_FAILURE);
               }
@@ -96,12 +95,11 @@ SteadyStateModel::checkPass(ModFileStructure &mod_file_struct, WarningConsolidat
     }
 
   set<int> orig_endogs = symbol_table.getOrigEndogenous();
-  for (set<int>::const_iterator it = orig_endogs.begin();
-       it != orig_endogs.end(); ++it)
+  for (int orig_endog : orig_endogs)
     {
-      if (find(so_far_defined.begin(), so_far_defined.end(), *it)
+      if (find(so_far_defined.begin(), so_far_defined.end(), orig_endog)
           == so_far_defined.end())
-        warnings << "WARNING: in the 'steady_state_model' block, variable '" << symbol_table.getName(*it) << "' is not assigned a value" << endl;
+        warnings << "WARNING: in the 'steady_state_model' block, variable '" << symbol_table.getName(orig_endog) << "' is not assigned a value" << endl;
     }
 }
 
@@ -135,12 +133,11 @@ SteadyStateModel::writeLatexSteadyStateFile(const string &basename) const
          << "\\begin{document}" << endl
          << "\\footnotesize" << endl;
 
-  for (vector<pair<vector<int>, expr_t> >::const_iterator it = def_table.begin();
-       it != def_table.end(); it++)
-    for (vector<int>::const_iterator it1 = it->first.begin(); it1 != it->first.end(); it1++)
+  for (const auto & it : def_table)
+    for (vector<int>::const_iterator it1 = it.first.begin(); it1 != it.first.end(); it1++)
       {
         int id = *it1;
-        expr_t value = it->second;
+        expr_t value = it.second;
         content_output << "\\begin{dmath}" << endl
                        << symbol_table.getTeXName(id) << " = ";
         value->writeOutput(content_output, oLatexStaticModel);
@@ -188,9 +185,9 @@ SteadyStateModel::writeSteadyStateFile(const string &basename, bool ramsey_model
            << "function steady_state!(ys_::Vector{Float64}, exo_::Vector{Float64}, "
            << "params::Vector{Float64})" << endl;
 
-  for (size_t i = 0; i < def_table.size(); i++)
+  for (const auto & i : def_table)
     {
-      const vector<int> &symb_ids = def_table[i].first;
+      const vector<int> &symb_ids = i.first;
       output << "    ";
       if (symb_ids.size() > 1)
         output << "[";
@@ -206,7 +203,7 @@ SteadyStateModel::writeSteadyStateFile(const string &basename, bool ramsey_model
         output << "]";
 
       output << "=";
-      def_table[i].second->writeOutput(output, output_type);
+      i.second->writeOutput(output, output_type);
       output << ";" << endl;
     }
   if (!julia)
@@ -251,9 +248,9 @@ SteadyStateModel::writeSteadyStateFileC(const string &basename, bool ramsey_mode
       return;
     }
 
-  for (size_t i = 0; i < def_table.size(); i++)
+  for (const auto & i : def_table)
     {
-      const vector<int> &symb_ids = def_table[i].first;
+      const vector<int> &symb_ids = i.first;
       output << "    ";
       if (symb_ids.size() > 1)
         std::cout << "Error: in C, multiple returns are not permitted in steady_state_model" << std::endl;
@@ -263,7 +260,7 @@ SteadyStateModel::writeSteadyStateFileC(const string &basename, bool ramsey_mode
         output << "double ";
       dynamic_cast<ExprNode *>(it->second)->writeOutput(output, oCSteadyStateFile);
       output << "=";
-      def_table[i].second->writeOutput(output, oCSteadyStateFile);
+      i.second->writeOutput(output, oCSteadyStateFile);
       output << ";" << endl;
     }
   output << "    // Auxiliary equations" << endl;
diff --git a/src/SymbolTable.cc b/src/SymbolTable.cc
index 396f8f88..c55498b3 100644
--- a/src/SymbolTable.cc
+++ b/src/SymbolTable.cc
@@ -71,10 +71,9 @@ SymbolTable::addSymbol(const string &name, SymbolType type, const string &tex_na
   string final_long_name = name;
   bool non_long_name_partition_exists = false;
   if (partition_value)
-    for (vector<pair<string *, string *> *>::const_iterator it = partition_value->begin();
-         it != partition_value->end(); it++)
-      if (*((*it)->first) == "long_name")
-        final_long_name = *((*it)->second);
+    for (auto it : *partition_value)
+      if (*(it->first) == "long_name")
+        final_long_name = *(it->second);
       else
         non_long_name_partition_exists = true;
 
@@ -88,9 +87,8 @@ SymbolTable::addSymbol(const string &name, SymbolType type, const string &tex_na
   if (non_long_name_partition_exists)
     {
       map<string, string> pmv;
-      for (vector<pair<string *, string *> *>::const_iterator it = partition_value->begin();
-           it != partition_value->end(); it++)
-        pmv[*((*it)->first)] = *((*it)->second);
+      for (auto it : *partition_value)
+        pmv[*(it->first)] = *(it->second);
       partition_value_map[id] = pmv;
     }
   return id;
@@ -198,15 +196,14 @@ map<string, map<int, string> >
 SymbolTable::getPartitionsForType(enum SymbolType st) const throw (UnknownSymbolIDException)
 {
   map<string, map<int, string> > partitions;
-  for (map<int, map<string, string> >::const_iterator it = partition_value_map.begin();
-       it != partition_value_map.end(); it++)
-    if (getType(it->first) == st)
-      for (map<string, string>::const_iterator it1 = it->second.begin();
-           it1 != it->second.end(); it1++)
+  for (const auto & it : partition_value_map)
+    if (getType(it.first) == st)
+      for (map<string, string>::const_iterator it1 = it.second.begin();
+           it1 != it.second.end(); it1++)
         {
           if (partitions.find(it1->first) == partitions.end())
             partitions[it1->first] = map<int, string> ();
-          partitions[it1->first][it->first] = it1->second;
+          partitions[it1->first][it.first] = it1->second;
         }
   return partitions;
 }
@@ -387,9 +384,8 @@ SymbolTable::writeOutput(ostream &output) const throw (NotYetFrozenException)
   if (predeterminedNbr() > 0)
     {
       output << "M_.predetermined_variables = [ ";
-      for (set<int>::const_iterator it = predetermined_variables.begin();
-           it != predetermined_variables.end(); it++)
-        output << getTypeSpecificID(*it)+1 << " ";
+      for (int predetermined_variable : predetermined_variables)
+        output << getTypeSpecificID(predetermined_variable)+1 << " ";
       output << "];" << endl;
     }
 
@@ -402,9 +398,8 @@ SymbolTable::writeOutput(ostream &output) const throw (NotYetFrozenException)
         output << "options_.varobs(" << ic << ")  = {'" << getName(*it) << "'};" << endl;
 
       output << "options_.varobs_id = [ ";
-      for (vector<int>::const_iterator it = varobs.begin();
-           it != varobs.end(); it++)
-        output << getTypeSpecificID(*it)+1 << " ";
+      for (int varob : varobs)
+        output << getTypeSpecificID(varob)+1 << " ";
       output << " ];"  << endl;
     }
 
@@ -417,9 +412,8 @@ SymbolTable::writeOutput(ostream &output) const throw (NotYetFrozenException)
         output << "options_.varexobs(" << ic << ")  = {'" << getName(*it) << "'};" << endl;
 
       output << "options_.varexobs_id = [ ";
-      for (vector<int>::const_iterator it = varexobs.begin();
-           it != varexobs.end(); it++)
-        output << getTypeSpecificID(*it)+1 << " ";
+      for (int varexob : varexobs)
+        output << getTypeSpecificID(varexob)+1 << " ";
       output << " ];"  << endl;
     }
 }
@@ -600,17 +594,14 @@ SymbolTable::writeCCOutput(ostream &output) const throw (NotYetFrozenException)
       output << "aux_vars.push_back(" << "av" << i << ");" << endl;
     }
 
-  for (set<int>::const_iterator it = predetermined_variables.begin();
-       it != predetermined_variables.end(); it++)
-    output << "predetermined_variables.push_back(" << getTypeSpecificID(*it) << ");" << endl;
+  for (int predetermined_variable : predetermined_variables)
+    output << "predetermined_variables.push_back(" << getTypeSpecificID(predetermined_variable) << ");" << endl;
 
-  for (vector<int>::const_iterator it = varobs.begin();
-       it != varobs.end(); it++)
-    output << "varobs.push_back(" << getTypeSpecificID(*it) << ");" << endl;
+  for (int varob : varobs)
+    output << "varobs.push_back(" << getTypeSpecificID(varob) << ");" << endl;
 
-  for (vector<int>::const_iterator it = varexobs.begin();
-       it != varexobs.end(); it++)
-    output << "varexobs.push_back(" << getTypeSpecificID(*it) << ");" << endl;
+  for (int varexob : varexobs)
+    output << "varexobs.push_back(" << getTypeSpecificID(varexob) << ");" << endl;
 }
 
 int
@@ -853,20 +844,20 @@ SymbolTable::addDiffForwardAuxiliaryVar(int orig_symb_id, expr_t expr_arg) throw
 int
 SymbolTable::searchAuxiliaryVars(int orig_symb_id, int orig_lead_lag) const throw (SearchFailedException)
 {
-  for (size_t i = 0; i < aux_vars.size(); i++)
-    if ((aux_vars[i].get_type() == avEndoLag || aux_vars[i].get_type() == avExoLag)
-        && aux_vars[i].get_orig_symb_id() == orig_symb_id && aux_vars[i].get_orig_lead_lag() == orig_lead_lag)
-      return aux_vars[i].get_symb_id();
+  for (const auto & aux_var : aux_vars)
+    if ((aux_var.get_type() == avEndoLag || aux_var.get_type() == avExoLag)
+        && aux_var.get_orig_symb_id() == orig_symb_id && aux_var.get_orig_lead_lag() == orig_lead_lag)
+      return aux_var.get_symb_id();
   throw SearchFailedException(orig_symb_id, orig_lead_lag);
 }
 
 int
 SymbolTable::getOrigSymbIdForAuxVar(int aux_var_symb_id) const throw (UnknownSymbolIDException)
 {
-  for (size_t i = 0; i < aux_vars.size(); i++)
-    if ((aux_vars[i].get_type() == avEndoLag || aux_vars[i].get_type() == avExoLag || aux_vars[i].get_type() == avDiff)
-        && aux_vars[i].get_symb_id() == aux_var_symb_id)
-      return aux_vars[i].get_orig_symb_id();
+  for (const auto & aux_var : aux_vars)
+    if ((aux_var.get_type() == avEndoLag || aux_var.get_type() == avExoLag || aux_var.get_type() == avDiff)
+        && aux_var.get_symb_id() == aux_var_symb_id)
+      return aux_var.get_orig_symb_id();
   throw UnknownSymbolIDException(aux_var_symb_id);
 }
 
@@ -874,10 +865,10 @@ expr_t
 SymbolTable::getAuxiliaryVarsExprNode(int symb_id) const throw (SearchFailedException)
 // throw exception if it is a Lagrange multiplier
 {
-  for (size_t i = 0; i < aux_vars.size(); i++)
-    if (aux_vars[i].get_symb_id() == symb_id)
+  for (const auto & aux_var : aux_vars)
+    if (aux_var.get_symb_id() == symb_id)
       {
-        expr_t expr_node = aux_vars[i].get_expr_node();
+        expr_t expr_node = aux_var.get_expr_node();
         if (expr_node != NULL)
           return expr_node;
         else
@@ -972,10 +963,9 @@ vector <int>
 SymbolTable::getTrendVarIds() const
 {
   vector <int> trendVars;
-  for (symbol_table_type::const_iterator it = symbol_table.begin();
-       it != symbol_table.end(); it++)
-    if (getType(it->second) == eTrend || getType(it->second) == eLogTrend)
-      trendVars.push_back(it->second);
+  for (const auto & it : symbol_table)
+    if (getType(it.second) == eTrend || getType(it.second) == eLogTrend)
+      trendVars.push_back(it.second);
   return trendVars;
 }
 
@@ -983,10 +973,9 @@ set<int>
 SymbolTable::getExogenous() const
 {
   set <int> exogs;
-  for (symbol_table_type::const_iterator it = symbol_table.begin();
-       it != symbol_table.end(); it++)
-    if (getType(it->second) == eExogenous)
-      exogs.insert(it->second);
+  for (const auto & it : symbol_table)
+    if (getType(it.second) == eExogenous)
+      exogs.insert(it.second);
   return exogs;
 }
 
@@ -994,11 +983,10 @@ set<int>
 SymbolTable::getObservedExogenous() const
 {
   set <int> oexogs;
-  for (symbol_table_type::const_iterator it = symbol_table.begin();
-       it != symbol_table.end(); it++)
-    if (getType(it->second) == eExogenous)
-      if (isObservedExogenousVariable(it->second))
-        oexogs.insert(it->second);
+  for (const auto & it : symbol_table)
+    if (getType(it.second) == eExogenous)
+      if (isObservedExogenousVariable(it.second))
+        oexogs.insert(it.second);
   return oexogs;
 }
 
@@ -1006,18 +994,17 @@ set<int>
 SymbolTable::getEndogenous() const
 {
   set <int> endogs;
-  for (symbol_table_type::const_iterator it = symbol_table.begin();
-       it != symbol_table.end(); it++)
-    if (getType(it->second) == eEndogenous)
-      endogs.insert(it->second);
+  for (const auto & it : symbol_table)
+    if (getType(it.second) == eEndogenous)
+      endogs.insert(it.second);
   return endogs;
 }
 
 bool
 SymbolTable::isAuxiliaryVariable(int symb_id) const
 {
-  for (int i = 0; i < (int) aux_vars.size(); i++)
-    if (aux_vars[i].get_symb_id() == symb_id)
+  for (const auto & aux_var : aux_vars)
+    if (aux_var.get_symb_id() == symb_id)
       return true;
   return false;
 }
@@ -1025,8 +1012,8 @@ SymbolTable::isAuxiliaryVariable(int symb_id) const
 bool
 SymbolTable::isAuxiliaryVariableButNotMultiplier(int symb_id) const
 {
-  for (int i = 0; i < (int) aux_vars.size(); i++)
-    if (aux_vars[i].get_symb_id() == symb_id && aux_vars[i].get_type() != avMultiplier)
+  for (const auto & aux_var : aux_vars)
+    if (aux_var.get_symb_id() == symb_id && aux_var.get_type() != avMultiplier)
       return true;
   return false;
 }
@@ -1035,10 +1022,9 @@ set<int>
 SymbolTable::getOrigEndogenous() const
 {
   set <int> origendogs;
-  for (symbol_table_type::const_iterator it = symbol_table.begin();
-       it != symbol_table.end(); it++)
-    if (getType(it->second) == eEndogenous && !isAuxiliaryVariable(it->second))
-      origendogs.insert(it->second);
+  for (const auto & it : symbol_table)
+    if (getType(it.second) == eEndogenous && !isAuxiliaryVariable(it.second))
+      origendogs.insert(it.second);
   return origendogs;
 }
 
@@ -1101,12 +1087,12 @@ SymbolTable::writeJuliaOutput(ostream &output) const throw (NotYetFrozenExceptio
     {
       output << "# Auxiliary Variables" << endl
              << "model_.aux_vars = [" << endl;
-      for (int i = 0; i < (int) aux_vars.size(); i++)
+      for (const auto & aux_var : aux_vars)
         {
           output << "                   DynareModel.AuxVars("
-                 << getTypeSpecificID(aux_vars[i].get_symb_id()) + 1 << ", "
-                 << aux_vars[i].get_type() << ", ";
-          switch (aux_vars[i].get_type())
+                 << getTypeSpecificID(aux_var.get_symb_id()) + 1 << ", "
+                 << aux_var.get_type() << ", ";
+          switch (aux_var.get_type())
             {
             case avEndoLead:
             case avExoLead:
@@ -1114,27 +1100,27 @@ SymbolTable::writeJuliaOutput(ostream &output) const throw (NotYetFrozenExceptio
             case avExoLag:
             case avVarModel:
             case avUnaryOp:
-              output << getTypeSpecificID(aux_vars[i].get_orig_symb_id()) + 1 << ", "
-                     << aux_vars[i].get_orig_lead_lag() << ", typemin(Int), string()";
+              output << getTypeSpecificID(aux_var.get_orig_symb_id()) + 1 << ", "
+                     << aux_var.get_orig_lead_lag() << ", typemin(Int), string()";
               break;
             case avDiff:
             case avDiffLag:
-              if (aux_vars[i].get_orig_symb_id() >= 0)
-                output << getTypeSpecificID(aux_vars[i].get_orig_symb_id()) + 1 << ", "
-                       << aux_vars[i].get_orig_lead_lag() << ", typemin(Int), string()";
+              if (aux_var.get_orig_symb_id() >= 0)
+                output << getTypeSpecificID(aux_var.get_orig_symb_id()) + 1 << ", "
+                       << aux_var.get_orig_lead_lag() << ", typemin(Int), string()";
               break;
             case avMultiplier:
-              output << "typemin(Int), typemin(Int), " << aux_vars[i].get_equation_number_for_multiplier() + 1
+              output << "typemin(Int), typemin(Int), " << aux_var.get_equation_number_for_multiplier() + 1
                      << ", string()";
               break;
             case avDiffForward:
-              output << getTypeSpecificID(aux_vars[i].get_orig_symb_id())+1 << ", typemin(Int), typemin(Int), string()";
+              output << getTypeSpecificID(aux_var.get_orig_symb_id())+1 << ", typemin(Int), typemin(Int), string()";
               break;
             case avExpectation:
               output << "typemin(Int), typemin(Int), typemin(Int), \"\\mathbb{E}_{t"
-                     << (aux_vars[i].get_information_set() < 0 ? "" : "+")
-                     << aux_vars[i].get_information_set() << "}(";
-              aux_vars[i].get_expr_node()->writeOutput(output, oLatexDynamicModel);
+                     << (aux_var.get_information_set() < 0 ? "" : "+")
+                     << aux_var.get_information_set() << "}(";
+              aux_var.get_expr_node()->writeOutput(output, oLatexDynamicModel);
               output << ")\"";
               break;
             default:
@@ -1149,10 +1135,9 @@ SymbolTable::writeJuliaOutput(ostream &output) const throw (NotYetFrozenExceptio
     {
       output << "# Predetermined Variables" << endl
              << "model_.pred_vars = [ " << endl;
-      for (set<int>::const_iterator it = predetermined_variables.begin();
-           it != predetermined_variables.end(); it++)
+      for (int predetermined_variable : predetermined_variables)
         output << "                   DynareModel.PredVars("
-               << getTypeSpecificID(*it)+1 << ")" << endl;
+               << getTypeSpecificID(predetermined_variable)+1 << ")" << endl;
       output << "                  ]" << endl;
     }
 
@@ -1160,10 +1145,9 @@ SymbolTable::writeJuliaOutput(ostream &output) const throw (NotYetFrozenExceptio
     {
       output << "# Observed Variables" << endl
              << "options_.obs_vars = [" << endl;
-      for (vector<int>::const_iterator it = varobs.begin();
-           it != varobs.end(); it++)
+      for (int varob : varobs)
         output << "                    DynareModel.ObsVars("
-               << getTypeSpecificID(*it)+1 << ")" << endl;
+               << getTypeSpecificID(varob)+1 << ")" << endl;
       output << "                   ]" << endl;
     }
 }
diff --git a/src/macro/MacroDriver.cc b/src/macro/MacroDriver.cc
index b0420b44..51543fe3 100644
--- a/src/macro/MacroDriver.cc
+++ b/src/macro/MacroDriver.cc
@@ -31,9 +31,8 @@ MacroDriver::MacroDriver()
 
 MacroDriver::~MacroDriver()
 {
-  for (set<const MacroValue *>::iterator it = values.begin();
-       it != values.end(); it++)
-    delete *it;
+  for (auto value : values)
+    delete value;
 }
 
 void
@@ -51,20 +50,19 @@ MacroDriver::parse(const string &f, const string &fb, const string &modfiletxt,
     an @#endif or an @#endfor - but no newline - no longer trigger an error.
   */
   stringstream file_with_endl;
-  for (map<string, string>::iterator it = defines.begin();
-       it != defines.end(); it++)
+  for (auto & define : defines)
     try
       {
-        boost::lexical_cast<int>(it->second);
-        file_with_endl << "@#define " << it->first << " = " << it->second << endl;
+        boost::lexical_cast<int>(define.second);
+        file_with_endl << "@#define " << define.first << " = " << define.second << endl;
       }
     catch (boost::bad_lexical_cast &)
       {
-        if (!it->second.empty() && it->second.at(0) == '[' && it->second.at(it->second.length()-1) == ']')
+        if (!define.second.empty() && define.second.at(0) == '[' && define.second.at(define.second.length()-1) == ']')
           // If the input is an array. Issue #1578
-          file_with_endl << "@#define " << it->first << " = " << it->second << endl;
+          file_with_endl << "@#define " << define.first << " = " << define.second << endl;
         else
-          file_with_endl << "@#define " << it->first << " = \"" << it->second << "\"" << endl;
+          file_with_endl << "@#define " << define.first << " = \"" << define.second << "\"" << endl;
       }
   file_with_endl << modfiletxt << endl;
 
@@ -220,9 +218,8 @@ MacroDriver::printvars(const Macro::parser::location_type &l, const bool tostdou
     {
       cout << "Macroprocessor: Printing macro variable values from " << file
            << " at line " << l.begin.line << endl;
-      for (map<string, const MacroValue *>::const_iterator it = env.begin();
-           it != env.end(); it++)
-        cout << "    " << it->first << " = " << it->second->print() << endl;
+      for (const auto & it : env)
+        cout << "    " << it.first << " = " << it.second->print() << endl;
       cout << endl;
       return "";
     }
@@ -231,8 +228,7 @@ MacroDriver::printvars(const Macro::parser::location_type &l, const bool tostdou
   if (!no_line_macro)
     intomfile << "@#line \"" << file << "\" " << l.begin.line << endl;
 
-  for (map<string, const MacroValue *>::const_iterator it = env.begin();
-       it != env.end(); it++)
-    intomfile<< "options_.macrovars_line_" << l.begin.line << "." << it->first << " = " << it->second->print() << ";" << endl;
+  for (const auto & it : env)
+    intomfile<< "options_.macrovars_line_" << l.begin.line << "." << it.first << " = " << it.second->print() << ";" << endl;
   return intomfile.str();
 }
diff --git a/src/macro/MacroValue.cc b/src/macro/MacroValue.cc
index db8c8ac2..950d276a 100644
--- a/src/macro/MacroValue.cc
+++ b/src/macro/MacroValue.cc
@@ -314,9 +314,8 @@ IntMV::in(const MacroValue *array) const throw (TypeError)
     throw TypeError("Type mismatch for 'in' operator");
 
   int result = 0;
-  for (vector<int>::const_iterator it = array2->values.begin();
-       it != array2->values.end(); it++)
-    if (*it == value)
+  for (int v : array2->values)
+    if (v == value)
       {
         result = 1;
         break;
@@ -387,12 +386,11 @@ StringMV::operator[](const MacroValue &mv) const throw (TypeError, OutOfBoundsEr
   if (mv2 == NULL)
     throw TypeError("Expression inside [] must be an integer array");
   string result;
-  for (vector<int>::const_iterator it = mv2->values.begin();
-       it != mv2->values.end(); it++)
+  for (int v : mv2->values)
     {
-      if (*it < 1 || *it > (int) value.length())
+      if (v < 1 || v > (int) value.length())
         throw OutOfBoundsError();
-      char c = value.at(*it - 1);
+      char c = value.at(v - 1);
       result.append(1, c);
     }
   return new StringMV(driver, result);
@@ -438,9 +436,8 @@ StringMV::in(const MacroValue *array) const throw (TypeError)
     throw TypeError("Type mismatch for 'in' operator");
 
   int result = 0;
-  for (vector<string>::const_iterator it = array2->values.begin();
-       it != array2->values.end(); it++)
-    if (*it == value)
+  for (const auto &v : array2->values)
+    if (v == value)
       {
         result = 1;
         break;
diff --git a/src/macro/MacroValue.hh b/src/macro/MacroValue.hh
index ced29189..ef1e8489 100644
--- a/src/macro/MacroValue.hh
+++ b/src/macro/MacroValue.hh
@@ -316,12 +316,11 @@ ArrayMV<T>::operator[](const MacroValue &mv) const throw (TypeError, OutOfBounds
   if (mv2 == NULL)
     throw TypeError("Expression inside [] must be an integer array");
   vector<T> result;
-  for (vector<int>::const_iterator it = mv2->values.begin();
-       it != mv2->values.end(); it++)
+  for (int value : mv2->values)
     {
-      if (*it < 1 || *it > (int) values.size())
+      if (value < 1 || value > (int) values.size())
         throw OutOfBoundsError();
-      result.push_back(values[*it - 1]);
+      result.push_back(values[value - 1]);
     }
 
   if (result.size() > 1 || result.size() == 0)