diff --git a/dynare++/integ/cc/quasi_mcarlo.cc b/dynare++/integ/cc/quasi_mcarlo.cc
index 131c398fe..2073c13af 100644
--- a/dynare++/integ/cc/quasi_mcarlo.cc
+++ b/dynare++/integ/cc/quasi_mcarlo.cc
@@ -122,8 +122,8 @@ HaltonSequence::operator=(const HaltonSequence &hs)
   num = hs.num;
   maxn = hs.maxn;
   ri.clear();
-  for (unsigned int i = 0; i < hs.ri.size(); i++)
-    ri.push_back(RadicalInverse(hs.ri[i]));
+  for (const auto & i : hs.ri)
+    ri.push_back(RadicalInverse(i));
   pt = hs.pt;
   return *this;
 }
@@ -134,8 +134,8 @@ HaltonSequence::operator=(const HaltonSequence &hs)
 void
 HaltonSequence::increase()
 {
-  for (unsigned int i = 0; i < ri.size(); i++)
-    ri[i].increase();
+  for (auto & i : ri)
+    i.increase();
   num++;
   if (num <= maxn)
     eval();
@@ -153,8 +153,8 @@ HaltonSequence::eval()
 void
 HaltonSequence::print() const
 {
-  for (unsigned int i = 0; i < ri.size(); i++)
-    ri[i].print();
+  for (const auto & i : ri)
+    i.print();
   printf("point=[ ");
   for (unsigned int i = 0; i < ri.size(); i++)
     printf("%7.6f ", pt[i]);
diff --git a/dynare++/integ/src/quadrature-points.cc b/dynare++/integ/src/quadrature-points.cc
index 56a4b7052..21604ae13 100644
--- a/dynare++/integ/src/quadrature-points.cc
+++ b/dynare++/integ/src/quadrature-points.cc
@@ -160,17 +160,17 @@ main(int argc, char **argv)
       // calculate weights and mass
       double mass = 0.0;
       std::vector<double> weights;
-      for (int i = 0; i < (int) points.size(); i++)
+      for (auto & point : points)
         {
-          weights.push_back(std::exp(-points[i]->dot(*(points[i]))));
+          weights.push_back(std::exp(-point->dot(*point)));
           mass += weights.back();
         }
 
       // calculate discarded mass
       double discard_mass = 0.0;
-      for (int i = 0; i < (int) weights.size(); i++)
-        if (weights[i]/mass < params.discard_weight)
-          discard_mass += weights[i];
+      for (double weight : weights)
+        if (weight/mass < params.discard_weight)
+          discard_mass += weight;
 
       printf("Total mass discarded:     %f\n", discard_mass/mass);
 
diff --git a/dynare++/kord/approximation.cc b/dynare++/kord/approximation.cc
index 82b85a2e3..c6b66e5a6 100644
--- a/dynare++/kord/approximation.cc
+++ b/dynare++/kord/approximation.cc
@@ -259,9 +259,9 @@ Approximation::saveRuleDerivs(const FGSContainer &g)
     }
   rule_ders = new FGSContainer(g);
   rule_ders_ss = new FGSContainer(4);
-  for (FGSContainer::iterator run = (*rule_ders).begin(); run != (*rule_ders).end(); ++run)
+  for (auto & run : (*rule_ders))
     {
-      FGSTensor *ten = new FGSTensor(ypart.nstat+ypart.npred, ypart.nyss(), *((*run).second));
+      FGSTensor *ten = new FGSTensor(ypart.nstat+ypart.npred, ypart.nyss(), *(run.second));
       rule_ders_ss->insert(ten);
     }
 }
diff --git a/dynare++/kord/decision_rule.cc b/dynare++/kord/decision_rule.cc
index 1287274fd..32d565ed2 100644
--- a/dynare++/kord/decision_rule.cc
+++ b/dynare++/kord/decision_rule.cc
@@ -32,10 +32,9 @@ FoldDecisionRule::FoldDecisionRule(const UnfoldDecisionRule &udr)
   : DecisionRuleImpl<KOrder::fold>(ctraits<KOrder::fold>::Tpol(udr.nrows(), udr.nvars()),
                                    udr.ypart, udr.nu, udr.ysteady)
 {
-  for (ctraits<KOrder::unfold>::Tpol::const_iterator it = udr.begin();
-       it != udr.end(); ++it)
+  for (const auto & it : udr)
     {
-      insert(new ctraits<KOrder::fold>::Ttensym(*((*it).second)));
+      insert(new ctraits<KOrder::fold>::Ttensym(*(it.second)));
     }
 }
 
@@ -44,10 +43,9 @@ UnfoldDecisionRule::UnfoldDecisionRule(const FoldDecisionRule &fdr)
   : DecisionRuleImpl<KOrder::unfold>(ctraits<KOrder::unfold>::Tpol(fdr.nrows(), fdr.nvars()),
                                      fdr.ypart, fdr.nu, fdr.ysteady)
 {
-  for (ctraits<KOrder::fold>::Tpol::const_iterator it = fdr.begin();
-       it != fdr.end(); ++it)
+  for (const auto & it : fdr)
     {
-      insert(new ctraits<KOrder::unfold>::Ttensym(*((*it).second)));
+      insert(new ctraits<KOrder::unfold>::Ttensym(*(it.second)));
     }
 }
 
@@ -196,11 +194,11 @@ SimResultsStats::calcMean()
   if (data.size()*num_per > 0)
     {
       double mult = 1.0/data.size()/num_per;
-      for (unsigned int i = 0; i < data.size(); i++)
+      for (auto & i : data)
         {
           for (int j = 0; j < num_per; j++)
             {
-              ConstVector col(*data[i], j);
+              ConstVector col(*i, j);
               mean.add(mult, col);
             }
         }
@@ -214,9 +212,9 @@ SimResultsStats::calcVcov()
     {
       vcov.zeros();
       double mult = 1.0/(data.size()*num_per - 1);
-      for (unsigned int i = 0; i < data.size(); i++)
+      for (auto & i : data)
         {
-          const TwoDMatrix &d = *(data[i]);
+          const TwoDMatrix &d = *i;
           for (int j = 0; j < num_per; j++)
             {
               for (int m = 0; m < num_y; m++)
@@ -276,9 +274,9 @@ SimResultsDynamicStats::calcMean()
       for (int j = 0; j < num_per; j++)
         {
           Vector meanj(mean, j);
-          for (unsigned int i = 0; i < data.size(); i++)
+          for (auto & i : data)
             {
-              ConstVector col(*data[i], j);
+              ConstVector col(*i, j);
               meanj.add(mult, col);
             }
         }
@@ -296,9 +294,9 @@ SimResultsDynamicStats::calcVariance()
         {
           ConstVector meanj(mean, j);
           Vector varj(variance, j);
-          for (int i = 0; i < (int) data.size(); i++)
+          for (auto & i : data)
             {
-              Vector col(ConstVector((*data[i]), j));
+              Vector col(ConstVector((*i), j));
               col.add(-1.0, meanj);
               for (int k = 0; k < col.length(); k++)
                 col[k] = col[k]*col[k];
@@ -350,8 +348,8 @@ SimResultsIRF::calcMeans()
   means.zeros();
   if (data.size() > 0)
     {
-      for (unsigned int i = 0; i < data.size(); i++)
-        means.add(1.0, *(data[i]));
+      for (auto & i : data)
+        means.add(1.0, *i);
       means.mult(1.0/data.size());
     }
 }
@@ -362,9 +360,9 @@ SimResultsIRF::calcVariances()
   if (data.size() > 1)
     {
       variances.zeros();
-      for (unsigned int i = 0; i < data.size(); i++)
+      for (auto & i : data)
         {
-          TwoDMatrix d((const TwoDMatrix &)(*(data[i])));
+          TwoDMatrix d((const TwoDMatrix &)(*i));
           d.add(-1.0, means);
           for (int j = 0; j < d.nrows(); j++)
             for (int k = 0; k < d.ncols(); k++)
@@ -450,9 +448,8 @@ IRFResults::IRFResults(const DynamicModel &mod, const DecisionRule &dr,
   pa << "Calculating IRFs against control for " << (int) irf_list_ind.size() << " shocks and for "
      << num_per << " periods" << endrec;
   const TwoDMatrix &vcov = mod.getVcov();
-  for (unsigned int ii = 0; ii < irf_list_ind.size(); ii++)
+  for (int ishock : irf_list_ind)
     {
-      int ishock = irf_list_ind[ii];
       double stderror = sqrt(vcov.get(ishock, ishock));
       irf_res.push_back(new SimResultsIRF(control, model.numeq(), num_per,
                                           ishock, stderror));
@@ -469,8 +466,8 @@ IRFResults::IRFResults(const DynamicModel &mod, const DecisionRule &dr,
 
 IRFResults::~IRFResults()
 {
-  for (unsigned int i = 0; i < irf_res.size(); i++)
-    delete irf_res[i];
+  for (auto & irf_re : irf_res)
+    delete irf_re;
 }
 
 void
diff --git a/dynare++/kord/journal.cc b/dynare++/kord/journal.cc
index 07231764e..395693afb 100644
--- a/dynare++/kord/journal.cc
+++ b/dynare++/kord/journal.cc
@@ -166,8 +166,8 @@ JournalRecord::operator<<(const IntSequence &s)
 void
 JournalRecord::writePrefix(const SystemResourcesFlash &f)
 {
-  for (int i = 0; i < MAXLEN; i++)
-    prefix[i] = ' ';
+  for (char & i : prefix)
+    i = ' ';
   double mb = 1024*1024;
   sprintf(prefix, "%07.6g", f.elapsed);
   sprintf(prefix+7, ":%c%05d", recChar, ord);
@@ -182,8 +182,8 @@ JournalRecord::writePrefix(const SystemResourcesFlash &f)
 void
 JournalRecordPair::writePrefixForEnd(const SystemResourcesFlash &f)
 {
-  for (int i = 0; i < MAXLEN; i++)
-    prefix_end[i] = ' ';
+  for (char & i : prefix_end)
+    i = ' ';
   double mb = 1024*1024;
   SystemResourcesFlash difnow;
   difnow.diff(f);
diff --git a/dynare++/parser/cc/atom_assignings.cc b/dynare++/parser/cc/atom_assignings.cc
index b1b51f65c..be1b919eb 100644
--- a/dynare++/parser/cc/atom_assignings.cc
+++ b/dynare++/parser/cc/atom_assignings.cc
@@ -18,9 +18,8 @@ AtomAssignings::AtomAssignings(const AtomAssignings &aa, ogp::StaticAtoms &a)
     order(aa.order)
 {
   // fill the lname2expr
-  for (Tvarintmap::const_iterator it = aa.lname2expr.begin();
-       it != aa.lname2expr.end(); ++it)
-    lname2expr.insert(Tvarintmap::value_type(left_names.query((*it).first), (*it).second));
+  for (auto it : aa.lname2expr)
+    lname2expr.insert(Tvarintmap::value_type(left_names.query(it.first), it.second));
 }
 
 /** A global symbol for passing info to the AtomAssignings from
@@ -153,11 +152,10 @@ AtomAssignings::apply_subst(const AtomSubstitutions::Toldnamemap &mm)
 {
   // go through all old variables and see what are their derived new
   // variables
-  for (AtomSubstitutions::Toldnamemap::const_iterator it = mm.begin();
-       it != mm.end(); ++it)
+  for (const auto & it : mm)
     {
-      const char *oldname = (*it).first;
-      const AtomSubstitutions::Tshiftnameset &sset = (*it).second;
+      const char *oldname = it.first;
+      const AtomSubstitutions::Tshiftnameset &sset = it.second;
       if (!sset.empty())
         {
           int told = atoms.index(oldname);
@@ -171,10 +169,9 @@ AtomAssignings::apply_subst(const AtomSubstitutions::Toldnamemap &mm)
           order.push_back(-1);
           // now go through all new names derived from the old name and
           // reference to the newly added formula
-          for (AtomSubstitutions::Tshiftnameset::const_iterator itt = sset.begin();
-               itt != sset.end(); ++itt)
+          for (const auto & itt : sset)
             {
-              const char *newname = (*itt).first;
+              const char *newname = itt.first;
               const char *nn = left_names.insert(newname);
               lname2expr.insert(Tvarintmap::value_type(nn, expr.nformulas()-1));
             }
@@ -188,9 +185,8 @@ AtomAssignings::print() const
   printf("Atom Assignings\nExpressions:\n");
   expr.print();
   printf("Left names:\n");
-  for (Tvarintmap::const_iterator it = lname2expr.begin();
-       it != lname2expr.end(); ++it)
-    printf("%s ==> %d (t=%d)\n", (*it).first, expr.formula((*it).second), order[(*it).second]);
+  for (auto it : lname2expr)
+    printf("%s ==> %d (t=%d)\n", it.first, expr.formula(it.second), order[it.second]);
 }
 
 void
diff --git a/dynare++/parser/cc/atom_substitutions.cc b/dynare++/parser/cc/atom_substitutions.cc
index 6c556ea71..c01bae638 100644
--- a/dynare++/parser/cc/atom_substitutions.cc
+++ b/dynare++/parser/cc/atom_substitutions.cc
@@ -14,20 +14,18 @@ AtomSubstitutions::AtomSubstitutions(const AtomSubstitutions &as, const FineAtom
   const NameStorage &ns = na.get_name_storage();
 
   // fill new2old
-  for (Tshiftmap::const_iterator it = as.new2old.begin();
-       it != as.new2old.end(); ++it)
-    new2old.insert(Tshiftmap::value_type(ns.query((*it).first),
-                                         Tshiftname(ns.query((*it).second.first),
-                                                    (*it).second.second)));
+  for (const auto & it : as.new2old)
+    new2old.insert(Tshiftmap::value_type(ns.query(it.first),
+                                         Tshiftname(ns.query(it.second.first),
+                                                    it.second.second)));
   // fill old2new
-  for (Toldnamemap::const_iterator it = as.old2new.begin();
-       it != as.old2new.end(); ++it)
+  for (const auto & it : as.old2new)
     {
       Tshiftnameset sset;
-      for (Tshiftnameset::const_iterator itt = (*it).second.begin();
-           itt != (*it).second.end(); ++itt)
+      for (Tshiftnameset::const_iterator itt = it.second.begin();
+           itt != it.second.end(); ++itt)
         sset.insert(Tshiftname(ns.query((*itt).first), (*itt).second));
-      old2new.insert(Toldnamemap::value_type(ns.query((*it).first), sset));
+      old2new.insert(Toldnamemap::value_type(ns.query(it.first), sset));
     }
 }
 
@@ -64,17 +62,15 @@ AtomSubstitutions::substitutions_finished(VarOrdering::ord_type ot)
   // create an external ordering of new_atoms from old_atoms
   const vector<const char *> &oa_ext = old_atoms.get_allvar();
   vector<const char *> na_ext;
-  for (unsigned int i = 0; i < oa_ext.size(); i++)
+  for (auto oname : oa_ext)
     {
-      const char *oname = oa_ext[i];
       // add the old name itself
       na_ext.push_back(oname);
       // add all new names derived from the old name
       Toldnamemap::const_iterator it = old2new.find(oname);
       if (it != old2new.end())
-        for (Tshiftnameset::const_iterator itt = (*it).second.begin();
-             itt != (*it).second.end(); ++itt)
-          na_ext.push_back((*itt).first);
+        for (const auto & itt : (*it).second)
+          na_ext.push_back(itt.first);
     }
 
   // call parsing finished for the new_atoms
@@ -88,10 +84,9 @@ AtomSubstitutions::get_new4old(const char *oldname, int tshift) const
   if (it != old2new.end())
     {
       const Tshiftnameset &sset = (*it).second;
-      for (Tshiftnameset::const_iterator itt = sset.begin();
-           itt != sset.end(); ++itt)
-        if ((*itt).second == -tshift)
-          return (*itt).first;
+      for (const auto & itt : sset)
+        if (itt.second == -tshift)
+          return itt.first;
     }
   return NULL;
 }
@@ -100,14 +95,14 @@ void
 AtomSubstitutions::print() const
 {
   printf("Atom Substitutions:\nOld ==> New:\n");
-  for (Toldnamemap::const_iterator it = old2new.begin(); it != old2new.end(); ++it)
-    for (Tshiftnameset::const_iterator itt = (*it).second.begin();
-         itt != (*it).second.end(); ++itt)
-      printf("    %s ==> [%s, %d]\n", (*it).first, (*itt).first, (*itt).second);
+  for (const auto & it : old2new)
+    for (Tshiftnameset::const_iterator itt = it.second.begin();
+         itt != it.second.end(); ++itt)
+      printf("    %s ==> [%s, %d]\n", it.first, (*itt).first, (*itt).second);
 
   printf("Old <== New:\n");
-  for (Tshiftmap::const_iterator it = new2old.begin(); it != new2old.end(); ++it)
-    printf("    [%s, %d] <== %s\n", (*it).second.first, (*it).second.second, (*it).first);
+  for (const auto & it : new2old)
+    printf("    [%s, %d] <== %s\n", it.second.first, it.second.second, it.first);
 }
 
 void
@@ -167,16 +162,14 @@ const char *
 SAtoms::findNameWithLeadInInterval(const vector<const char *> &names,
                                    int ll1, int ll2) const
 {
-  for (unsigned int i = 0; i < names.size(); i++)
+  for (auto name : names)
     {
-      const char *name = names[i];
       DynamicAtoms::Tvarmap::const_iterator it = vars.find(name);
       if (it != vars.end())
         {
           const DynamicAtoms::Tlagmap &lmap = (*it).second;
-          for (DynamicAtoms::Tlagmap::const_iterator itt = lmap.begin();
-               itt != lmap.end(); ++itt)
-            if ((*itt).first >= ll1 && (*itt).first <= ll2)
+          for (auto itt : lmap)
+            if (itt.first >= ll1 && itt.first <= ll2)
               return name;
         }
     }
diff --git a/dynare++/parser/cc/dynamic_atoms.cc b/dynare++/parser/cc/dynamic_atoms.cc
index d8f424109..77f78e200 100644
--- a/dynare++/parser/cc/dynamic_atoms.cc
+++ b/dynare++/parser/cc/dynamic_atoms.cc
@@ -12,10 +12,10 @@ using namespace ogp;
 
 NameStorage::NameStorage(const NameStorage &stor)
 {
-  for (unsigned int i = 0; i < stor.name_store.size(); i++)
+  for (auto i : stor.name_store)
     {
-      char *str = new char[strlen(stor.name_store[i])+1];
-      strcpy(str, stor.name_store[i]);
+      char *str = new char[strlen(i)+1];
+      strcpy(str, i);
       name_store.push_back(str);
       name_set.insert(str);
     }
@@ -61,20 +61,19 @@ NameStorage::insert(const char *name)
 void
 NameStorage::print() const
 {
-  for (unsigned int i = 0; i < name_store.size(); i++)
-    printf("%s\n", name_store[i]);
+  for (auto i : name_store)
+    printf("%s\n", i);
 }
 
 void
 Constants::import_constants(const Constants &c, OperationTree &otree, Tintintmap &tmap)
 {
-  for (Tconstantmap::const_iterator it = c.cmap.begin();
-       it != c.cmap.end(); ++it)
+  for (auto it : c.cmap)
     {
-      int told = (*it).first;
+      int told = it.first;
       int tnew = otree.add_nulary();
       tmap.insert(Tintintmap::value_type(told, tnew));
-      add_constant(tnew, (*it).second);
+      add_constant(tnew, it.second);
     }
 }
 
@@ -147,15 +146,13 @@ DynamicAtoms::DynamicAtoms(const DynamicAtoms &da)
     nv(da.nv), minlag(da.minlag), maxlead(da.maxlead)
 {
   // copy vars
-  for (Tvarmap::const_iterator it = da.vars.begin();
-       it != da.vars.end(); ++it)
-    vars.insert(Tvarmap::value_type(varnames.query((*it).first),
-                                    (*it).second));
+  for (const auto & var : da.vars)
+    vars.insert(Tvarmap::value_type(varnames.query(var.first),
+                                    var.second));
   // copy indices
-  for (Tindexmap::const_iterator it = da.indices.begin();
-       it != da.indices.end(); ++it)
-    indices.insert(Tindexmap::value_type((*it).first,
-                                         varnames.query((*it).second)));
+  for (auto indice : da.indices)
+    indices.insert(Tindexmap::value_type(indice.first,
+                                         varnames.query(indice.second)));
 }
 
 int
@@ -294,9 +291,9 @@ DynamicAtoms::update_minmaxll()
   for (Tvarmap::const_iterator it = vars.begin(); it != vars.end(); ++it)
     {
       const Tlagmap &lmap = (*it).second;
-      for (Tlagmap::const_iterator itt = lmap.begin(); itt != lmap.end(); ++itt)
+      for (auto itt : lmap)
         {
-          int ll = (*itt).first;
+          int ll = itt.first;
           if (ll < minlag)
             minlag = ll;
           if (ll > maxlead)
@@ -309,13 +306,11 @@ vector<int>
 DynamicAtoms::variables() const
 {
   vector<int> res;
-  for (Tvarmap::const_iterator it = vars.begin();
-       it != vars.end(); ++it)
+  for (const auto & var : vars)
     {
-      const Tlagmap &lmap = (*it).second;
-      for (Tlagmap::const_iterator itt = lmap.begin();
-           itt != lmap.end(); ++itt)
-        res.push_back((*itt).second);
+      const Tlagmap &lmap = var.second;
+      for (auto itt : lmap)
+        res.push_back(itt.second);
     }
   return res;
 }
@@ -355,10 +350,10 @@ DynamicAtoms::varspan(const vector<const char *> &names, int &mlead, int &mlag)
 {
   mlead = INT_MIN;
   mlag = INT_MAX;
-  for (unsigned int i = 0; i < names.size(); i++)
+  for (auto name : names)
     {
       int lag, lead;
-      varspan(names[i], lead, lag);
+      varspan(name, lead, lag);
       if (lead > mlead)
         mlead = lead;
       if (lag < mlag)
@@ -436,18 +431,15 @@ DynamicAtoms::print() const
   printf("constants:\n");
   Constants::print();
   printf("variables:\n");
-  for (Tvarmap::const_iterator it = vars.begin();
-       it != vars.end(); ++it)
+  for (const auto & var : vars)
     {
-      const Tlagmap &lmap = (*it).second;
-      for (Tlagmap::const_iterator itt = lmap.begin();
-           itt != lmap.end(); ++itt)
-        printf("$%d: %s(%d)\n", (*itt).second, (*it).first, (*itt).first);
+      const Tlagmap &lmap = var.second;
+      for (auto itt : lmap)
+        printf("$%d: %s(%d)\n", itt.second, var.first, itt.first);
     }
   printf("indices:\n");
-  for (Tindexmap::const_iterator it = indices.begin();
-       it != indices.end(); ++it)
-    printf("t=%d ==> %s\n", (*it).first, (*it).second);
+  for (auto indice : indices)
+    printf("t=%d ==> %s\n", indice.first, indice.second);
 }
 
 /** Note that the str has been parsed by the lexicographic
@@ -578,12 +570,12 @@ VarOrdering::do_general(ord_type ordering)
 
   // make der_atoms and positions
   int off = 0;
-  for (unsigned int i = 0; i < 8; i++)
-    for (unsigned int j = 0; j < (ords[i])->size(); j++, off++)
-      if ((*(ords[i]))[j] != -1)
+  for (auto & ord : ords)
+    for (unsigned int j = 0; j < ord->size(); j++, off++)
+      if ((*ord)[j] != -1)
         {
-          der_atoms.push_back((*(ords[i]))[j]);
-          positions.insert(std::pair<int, int>((*(ords[i]))[j], off));
+          der_atoms.push_back((*ord)[j]);
+          positions.insert(std::pair<int, int>((*ord)[j], off));
         }
 
   // set integer constants
@@ -618,11 +610,10 @@ VarOrdering::do_increasing_time()
       try
         {
           const DynamicAtoms::Tlagmap &lmap = atoms.lagmap(varnames[iv]);
-          for (DynamicAtoms::Tlagmap::const_iterator it = lmap.begin();
-               it != lmap.end(); ++it)
+          for (auto it : lmap)
             {
-              int ll = (*it).first;
-              int t = (*it).second;
+              int ll = it.first;
+              int t = it.second;
               tree_ind[ll-mlag][iv] = t;
             }
         }
@@ -653,10 +644,10 @@ VarOrdering::do_increasing_time()
     }
 
   // set n_stat, n_pred, n_both, and n_forw
-  for (unsigned int iv = 0; iv < varnames.size(); iv++)
+  for (auto varname : varnames)
     {
       int mmlag, mmlead;
-      atoms.varspan(varnames[iv], mmlead, mmlag);
+      atoms.varspan(varname, mmlead, mmlag);
       if (mmlead == 0 && mmlag == 0)
         {
           n_stat++;
@@ -691,17 +682,17 @@ VarOrdering::print() const
 {
   printf("nstat=%d, npred=%d, nboth=%d, nforw=%d\n", n_stat, n_pred, n_both, n_forw);
   printf("der_atoms:\n");
-  for (unsigned int i = 0; i < der_atoms.size(); i++)
-    printf(" %d", der_atoms[i]);
+  for (int der_atom : der_atoms)
+    printf(" %d", der_atom);
   printf("\nmap:\n");
-  for (map<int, int>::const_iterator it = positions.begin(); it != positions.end(); ++it)
-    printf(" [%d->%d]", (*it).first, (*it).second);
+  for (auto position : positions)
+    printf(" [%d->%d]", position.first, position.second);
   printf("\ny2outer:\n");
-  for (unsigned int i = 0; i < y2outer.size(); i++)
-    printf(" %d", y2outer[i]);
+  for (int i : y2outer)
+    printf(" %d", i);
   printf("\nouter2y:\n");
-  for (unsigned int i = 0; i < outer2y.size(); i++)
-    printf(" %d", outer2y[i]);
+  for (int i : outer2y)
+    printf(" %d", i);
   printf("\n");
 }
 
diff --git a/dynare++/parser/cc/fine_atoms.cc b/dynare++/parser/cc/fine_atoms.cc
index 2080452f8..3b747653d 100644
--- a/dynare++/parser/cc/fine_atoms.cc
+++ b/dynare++/parser/cc/fine_atoms.cc
@@ -16,14 +16,14 @@ AllvarOuterOrdering::AllvarOuterOrdering(const vector<const char *> &allvar_oute
     exo2all(a.get_exovars().size(), -1)
 {
   // fill in the allvar from allvar_outer
-  for (unsigned int i = 0; i < allvar_outer.size(); i++)
+  for (auto i : allvar_outer)
     {
-      const char *s = atoms.varnames.query(allvar_outer[i]);
+      const char *s = atoms.varnames.query(i);
       if (s)
         allvar.push_back(s);
       else
         throw ogu::Exception(__FILE__, __LINE__,
-                             string("Variable ") + allvar_outer[i] + " is not a declared symbol in AllvarOuterOrdering constructor");
+                             string("Variable ") + i + " is not a declared symbol in AllvarOuterOrdering constructor");
     }
 
   // fill in endo2all and exo2all
@@ -67,9 +67,9 @@ AllvarOuterOrdering::AllvarOuterOrdering(const AllvarOuterOrdering &avo,
     exo2all(avo.exo2all)
 {
   // fill in the allvar from avo.allvar
-  for (unsigned int i = 0; i < avo.allvar.size(); i++)
+  for (auto i : avo.allvar)
     {
-      const char *s = atoms.varnames.query(avo.allvar[i]);
+      const char *s = atoms.varnames.query(i);
       allvar.push_back(s);
     }
 }
@@ -82,32 +82,32 @@ FineAtoms::FineAtoms(const FineAtoms &fa)
     exo_atoms_map(fa.exo_atoms_map)
 {
   // fill in params
-  for (unsigned int i = 0; i < fa.params.size(); i++)
+  for (auto param : fa.params)
     {
-      const char *s = varnames.query(fa.params[i]);
+      const char *s = varnames.query(param);
       if (!s)
         throw ogu::Exception(__FILE__, __LINE__,
-                             string("Parameter ") + fa.params[i] + " does not exist in FineAtoms copy cosntructor");
+                             string("Parameter ") + param + " does not exist in FineAtoms copy cosntructor");
       params.push_back(s);
       param_outer_map.insert(Tvarintmap::value_type(s, params.size()-1));
     }
   // fill in endovars
-  for (unsigned int i = 0; i < fa.endovars.size(); i++)
+  for (auto endovar : fa.endovars)
     {
-      const char *s = varnames.query(fa.endovars[i]);
+      const char *s = varnames.query(endovar);
       if (!s)
         throw ogu::Exception(__FILE__, __LINE__,
-                             string("Endo variable ") + fa.endovars[i] + " does not exist in FineAtoms copy constructor");
+                             string("Endo variable ") + endovar + " does not exist in FineAtoms copy constructor");
       endovars.push_back(s);
       endo_outer_map.insert(Tvarintmap::value_type(s, endovars.size()-1));
     }
   // fill in exovars
-  for (unsigned int i = 0; i < fa.exovars.size(); i++)
+  for (auto exovar : fa.exovars)
     {
-      const char *s = varnames.query(fa.exovars[i]);
+      const char *s = varnames.query(exovar);
       if (!s)
         throw ogu::Exception(__FILE__, __LINE__,
-                             string("Exo variable ") + fa.exovars[i] + " does not exist in FineAtoms copy cosntructor");
+                             string("Exo variable ") + exovar + " does not exist in FineAtoms copy cosntructor");
       exovars.push_back(s);
       exo_outer_map.insert(Tvarintmap::value_type(s, exovars.size()-1));
     }
diff --git a/dynare++/parser/cc/formula_parser.cc b/dynare++/parser/cc/formula_parser.cc
index 58e674390..fefd1e424 100644
--- a/dynare++/parser/cc/formula_parser.cc
+++ b/dynare++/parser/cc/formula_parser.cc
@@ -20,8 +20,8 @@ FormulaParser::FormulaParser(const FormulaParser &fp, Atoms &a)
   : otree(fp.otree), atoms(a), formulas(fp.formulas), ders()
 {
   // create derivatives
-  for (unsigned int i = 0; i < fp.ders.size(); i++)
-    ders.push_back(new FormulaDerivatives(*(fp.ders[i])));
+  for (auto der : fp.ders)
+    ders.push_back(new FormulaDerivatives(*der));
 }
 
 FormulaParser::~FormulaParser()
@@ -35,8 +35,8 @@ FormulaParser::differentiate(int max_order)
   destroy_derivatives();
   vector<int> vars;
   vars = atoms.variables();
-  for (unsigned int i = 0; i < formulas.size(); i++)
-    ders.push_back(new FormulaDerivatives(otree, vars, formulas[i], max_order));
+  for (int formula : formulas)
+    ders.push_back(new FormulaDerivatives(otree, vars, formula, max_order));
 }
 
 const FormulaDerivatives &
@@ -156,9 +156,9 @@ int
 FormulaParser::last_formula() const
 {
   int res = -1;
-  for (unsigned int i = 0; i < formulas.size(); i++)
-    if (res < formulas[i])
-      res = formulas[i];
+  for (int formula : formulas)
+    if (res < formula)
+      res = formula;
   return std::max(res, otree.get_last_nulary());
 }
 
@@ -181,10 +181,10 @@ void
 FormulaParser::print() const
 {
   atoms.print();
-  for (unsigned int i = 0; i < formulas.size(); i++)
+  for (int formula : formulas)
     {
-      printf("formula %d:\n", formulas[i]);
-      otree.print_operation(formulas[i]);
+      printf("formula %d:\n", formula);
+      otree.print_operation(formula);
     }
   for (unsigned int i = 0; i < ders.size(); i++)
     {
@@ -283,13 +283,12 @@ FormulaDerivatives::derivative(const FoldMultiIndex &mi) const
 void
 FormulaDerivatives::print(const OperationTree &otree) const
 {
-  for (Tfmiintmap::const_iterator it = ind2der.begin();
-       it != ind2der.end(); ++it)
+  for (const auto & it : ind2der)
     {
       printf("derivative ");
-      (*it).first.print();
-      printf(" is formula %d\n", tder[(*it).second]);
-      otree.print_operation(tder[(*it).second]);
+      it.first.print();
+      printf(" is formula %d\n", tder[it.second]);
+      otree.print_operation(tder[it.second]);
     }
 }
 
@@ -480,8 +479,8 @@ ltfmi::operator()(const FoldMultiIndex &i1, const FoldMultiIndex &i2) const
 FormulaDerEvaluator::FormulaDerEvaluator(const FormulaParser &fp)
   : etree(fp.otree, -1)
 {
-  for (unsigned int i = 0; i < fp.ders.size(); i++)
-    ders.push_back((const FormulaDerivatives *) (fp.ders[i]));
+  for (auto der : fp.ders)
+    ders.push_back((const FormulaDerivatives *) der);
 
   der_atoms = fp.atoms.variables();
 }
diff --git a/dynare++/parser/cc/static_atoms.cc b/dynare++/parser/cc/static_atoms.cc
index 288f2a37a..2a045ccd2 100644
--- a/dynare++/parser/cc/static_atoms.cc
+++ b/dynare++/parser/cc/static_atoms.cc
@@ -12,26 +12,24 @@ StaticAtoms::StaticAtoms(const StaticAtoms &a)
     varorder(), vars(), indices()
 {
   // fill varorder
-  for (unsigned int i = 0; i < a.varorder.size(); i++)
+  for (auto i : a.varorder)
     {
-      const char *s = varnames.query(a.varorder[i]);
+      const char *s = varnames.query(i);
       varorder.push_back(s);
     }
 
   // fill vars
-  for (Tvarmap::const_iterator it = a.vars.begin();
-       it != a.vars.end(); ++it)
+  for (auto var : a.vars)
     {
-      const char *s = varnames.query((*it).first);
-      vars.insert(Tvarmap::value_type(s, (*it).second));
+      const char *s = varnames.query(var.first);
+      vars.insert(Tvarmap::value_type(s, var.second));
     }
 
   // fill indices
-  for (Tinvmap::const_iterator it = a.indices.begin();
-       it != a.indices.end(); ++it)
+  for (auto indice : a.indices)
     {
-      const char *s = varnames.query((*it).second);
-      indices.insert(Tinvmap::value_type((*it).first, s));
+      const char *s = varnames.query(indice.second);
+      indices.insert(Tinvmap::value_type(indice.first, s));
     }
 }
 
@@ -49,10 +47,9 @@ StaticAtoms::import_atoms(const DynamicAtoms &da, OperationTree &otree, Tintintm
       if (da.is_referenced(name))
         {
           const DynamicAtoms::Tlagmap &lmap = da.lagmap(name);
-          for (DynamicAtoms::Tlagmap::const_iterator it = lmap.begin();
-               it != lmap.end(); ++it)
+          for (auto it : lmap)
             {
-              int told = (*it).second;
+              int told = it.second;
               tmap.insert(Tintintmap::value_type(told, tnew));
             }
         }
@@ -113,10 +110,9 @@ vector<int>
 StaticAtoms::variables() const
 {
   vector<int> res;
-  for (Tvarmap::const_iterator it = vars.begin();
-       it != vars.end(); ++it)
+  for (auto var : vars)
     {
-      res.push_back((*it).second);
+      res.push_back(var.second);
     }
   return res;
 }
@@ -136,6 +132,6 @@ StaticAtoms::print() const
   printf("variable names:\n");
   varnames.print();
   printf("map to tree indices:\n");
-  for (Tvarmap::const_iterator it = vars.begin(); it != vars.end(); ++it)
-    printf("%s\t->\t%d\n", (*it).first, (*it).second);
+  for (auto var : vars)
+    printf("%s\t->\t%d\n", var.first, var.second);
 }
diff --git a/dynare++/parser/cc/static_fine_atoms.cc b/dynare++/parser/cc/static_fine_atoms.cc
index 807e9df8c..8b37e3215 100644
--- a/dynare++/parser/cc/static_fine_atoms.cc
+++ b/dynare++/parser/cc/static_fine_atoms.cc
@@ -50,18 +50,18 @@ StaticFineAtoms::import_atoms(const FineAtoms &fa, OperationTree &otree, Tintint
 
   // parameters
   const vector<const char *> &fa_params = fa.get_params();
-  for (unsigned int i = 0; i < fa_params.size(); i++)
-    register_param(fa_params[i]);
+  for (auto fa_param : fa_params)
+    register_param(fa_param);
 
   // endogenous
   const vector<const char *> &fa_endovars = fa.get_endovars();
-  for (unsigned int i = 0; i < fa_endovars.size(); i++)
-    register_endo(fa_endovars[i]);
+  for (auto fa_endovar : fa_endovars)
+    register_endo(fa_endovar);
 
   // exogenous
   const vector<const char *> &fa_exovars = fa.get_exovars();
-  for (unsigned int i = 0; i < fa_exovars.size(); i++)
-    register_exo(fa_exovars[i]);
+  for (auto fa_exovar : fa_exovars)
+    register_exo(fa_exovar);
 
   parsing_finished();
 }
@@ -77,8 +77,8 @@ StaticFineAtoms::import_atoms(const FineAtoms &fa, OperationTree &otree, Tintint
 
   // parameters
   const vector<const char *> &fa_params = fa.get_params();
-  for (unsigned int i = 0; i < fa_params.size(); i++)
-    register_param(fa_params[i]);
+  for (auto fa_param : fa_params)
+    register_param(fa_param);
 
   // endogenous
   const vector<const char *> &fa_endovars = fa.get_endovars();
@@ -112,18 +112,18 @@ StaticFineAtoms::parsing_finished()
 
   // go through all endo and exo insert tree indices, ignore names
   // whose tree index is -1 (those which are not referenced)
-  for (unsigned int i = 0; i < endovars.size(); i++)
+  for (auto & endovar : endovars)
     {
-      int t = index(endovars[i]);
+      int t = index(endovar);
       if (t != -1)
         {
           endo_atoms_map.push_back(der_atoms.size());
           der_atoms.push_back(t);
         }
     }
-  for (unsigned int i = 0; i < exovars.size(); i++)
+  for (auto & exovar : exovars)
     {
-      int t = index(exovars[i]);
+      int t = index(exovar);
       if (t != -1)
         {
           exo_atoms_map.push_back(der_atoms.size());
diff --git a/dynare++/parser/cc/tree.cc b/dynare++/parser/cc/tree.cc
index 8209e9c0c..a6bc9dbd9 100644
--- a/dynare++/parser/cc/tree.cc
+++ b/dynare++/parser/cc/tree.cc
@@ -468,8 +468,8 @@ OperationTree::select_terms_inv(int t, const opselector &sel, unordered_set<int>
 void
 OperationTree::forget_derivative_maps()
 {
-  for (unsigned int i = 0; i < derivatives.size(); i++)
-    derivatives[i].clear();
+  for (auto & derivative : derivatives)
+    derivative.clear();
 }
 
 void
diff --git a/dynare++/src/dynare3.cc b/dynare++/src/dynare3.cc
index 80826598e..08603cab2 100644
--- a/dynare++/src/dynare3.cc
+++ b/dynare++/src/dynare3.cc
@@ -21,14 +21,14 @@ vector<int>
 DynareNameList::selectIndices(const vector<const char *> &ns) const
 {
   vector<int> res;
-  for (unsigned int i = 0; i < ns.size(); i++)
+  for (auto n : ns)
     {
       int j = 0;
-      while (j < getNum() && strcmp(getName(j), ns[i]) != 0)
+      while (j < getNum() && strcmp(getName(j), n) != 0)
         j++;
       if (j == getNum())
         throw DynareException(__FILE__, __LINE__,
-                              string("Couldn't find name for ") + ns[i]
+                              string("Couldn't find name for ") + n
                               +" in DynareNameList::selectIndices");
       res.push_back(j);
     }
diff --git a/dynare++/src/dynare_atoms.cc b/dynare++/src/dynare_atoms.cc
index 026631d9f..68b952eda 100644
--- a/dynare++/src/dynare_atoms.cc
+++ b/dynare++/src/dynare_atoms.cc
@@ -37,9 +37,8 @@ DynareDynamicAtoms::DynareDynamicAtoms(const DynareDynamicAtoms &dda)
   : SAtoms(dda)
 {
   // fill atom_type
-  for (Tatypemap::const_iterator it = dda.atom_type.begin();
-       it != dda.atom_type.end(); ++it)
-    atom_type.insert(Tatypemap::value_type(varnames.query((*it).first), (*it).second));
+  for (auto it : dda.atom_type)
+    atom_type.insert(Tatypemap::value_type(varnames.query(it.first), it.second));
 }
 
 void
@@ -101,10 +100,9 @@ DynareDynamicAtoms::print() const
 {
   SAtoms::print();
   printf("Name types:\n");
-  for (Tatypemap::const_iterator it = atom_type.begin();
-       it != atom_type.end(); ++it)
-    printf("name=%s type=%s\n", (*it).first,
-           ((*it).second == endovar) ? "endovar" : (((*it).second == exovar) ? "exovar" : "param"));
+  for (auto it : atom_type)
+    printf("name=%s type=%s\n", it.first,
+           (it.second == endovar) ? "endovar" : ((it.second == exovar) ? "exovar" : "param"));
 }
 
 std::string
@@ -154,10 +152,9 @@ DynareAtomValues::setValues(ogp::EvalTree &et) const
       if (atoms.is_referenced(atoms.get_params()[i]))
         {
           const ogp::DynamicAtoms::Tlagmap &lmap = atoms.lagmap(atoms.get_params()[i]);
-          for (ogp::DynamicAtoms::Tlagmap::const_iterator it = lmap.begin();
-               it != lmap.end(); ++it)
+          for (auto it : lmap)
             {
-              int t = (*it).second;
+              int t = it.second;
               et.set_nulary(t, paramvals[i]);
             }
         }
@@ -169,11 +166,10 @@ DynareAtomValues::setValues(ogp::EvalTree &et) const
       if (atoms.is_referenced(atoms.get_endovars()[outer_i]))
         {
           const ogp::DynamicAtoms::Tlagmap &lmap = atoms.lagmap(atoms.get_endovars()[outer_i]);
-          for (ogp::DynamicAtoms::Tlagmap::const_iterator it = lmap.begin();
-               it != lmap.end(); ++it)
+          for (auto it : lmap)
             {
-              int ll = (*it).first;
-              int t = (*it).second;
+              int ll = it.first;
+              int t = it.second;
               int i = atoms.outer2y_endo()[outer_i];
               if (ll == -1)
                 {
@@ -193,13 +189,12 @@ DynareAtomValues::setValues(ogp::EvalTree &et) const
       if (atoms.is_referenced(atoms.get_exovars()[outer_i]))
         {
           const ogp::DynamicAtoms::Tlagmap &lmap = atoms.lagmap(atoms.get_exovars()[outer_i]);
-          for (ogp::DynamicAtoms::Tlagmap::const_iterator it = lmap.begin();
-               it != lmap.end(); ++it)
+          for (auto it : lmap)
             {
-              int ll = (*it).first;
+              int ll = it.first;
               if (ll == 0)   // this is always true because of checks
                 {
-                  int t = (*it).second;
+                  int t = it.second;
                   int i = atoms.outer2y_exo()[outer_i];
                   et.set_nulary(t, xx[i]);
                 }
@@ -215,9 +210,8 @@ DynareStaticSteadyAtomValues::setValues(ogp::EvalTree &et) const
   atoms_static.setValues(et);
 
   // set parameters
-  for (unsigned int i = 0; i < atoms_static.get_params().size(); i++)
+  for (auto name : atoms_static.get_params())
     {
-      const char *name = atoms_static.get_params()[i];
       int t = atoms_static.index(name);
       if (t != -1)
         {
@@ -227,9 +221,8 @@ DynareStaticSteadyAtomValues::setValues(ogp::EvalTree &et) const
     }
 
   // set endogenous
-  for (unsigned int i = 0; i < atoms_static.get_endovars().size(); i++)
+  for (auto name : atoms_static.get_endovars())
     {
-      const char *name = atoms_static.get_endovars()[i];
       int t = atoms_static.index(name);
       if (t != -1)
         {
@@ -239,9 +232,8 @@ DynareStaticSteadyAtomValues::setValues(ogp::EvalTree &et) const
     }
 
   // set exogenous
-  for (unsigned int i = 0; i < atoms_static.get_exovars().size(); i++)
+  for (auto name : atoms_static.get_exovars())
     {
-      const char *name = atoms_static.get_exovars()[i];
       int t = atoms_static.index(name);
       if (t != -1)
         et.set_nulary(t, 0.0);
@@ -255,11 +247,10 @@ DynareSteadySubstitutions::DynareSteadySubstitutions(const ogp::FineAtoms &a,
   : atoms(a), y(yy)
 {
   // fill the vector of left and right hand sides
-  for (Tsubstmap::const_iterator it = subst.begin();
-       it != subst.end(); ++it)
+  for (auto it : subst)
     {
-      left_hand_sides.push_back((*it).first);
-      right_hand_sides.push_back((*it).second);
+      left_hand_sides.push_back(it.first);
+      right_hand_sides.push_back(it.second);
     }
 
   // evaluate right hand sides
@@ -286,11 +277,10 @@ DynareStaticSteadySubstitutions(const ogp::FineAtoms &a, const ogp::StaticFineAt
   : atoms(a), atoms_static(sa), y(yy)
 {
   // fill the vector of left and right hand sides
-  for (Tsubstmap::const_iterator it = subst.begin();
-       it != subst.end(); ++it)
+  for (auto it : subst)
     {
-      left_hand_sides.push_back((*it).first);
-      right_hand_sides.push_back((*it).second);
+      left_hand_sides.push_back(it.first);
+      right_hand_sides.push_back(it.second);
     }
 
   // evaluate right hand sides
diff --git a/dynare++/src/dynare_model.cc b/dynare++/src/dynare_model.cc
index 7038f118c..89ba8b4d5 100644
--- a/dynare++/src/dynare_model.cc
+++ b/dynare++/src/dynare_model.cc
@@ -130,20 +130,20 @@ DynareModel::dump_model(std::ostream &os) const
 {
   // endogenous variable declaration
   os << "var";
-  for (int i = 0; i < (int) atoms.get_endovars().size(); i++)
-    os << " " << atoms.get_endovars()[i];
+  for (auto i : atoms.get_endovars())
+    os << " " << i;
   os << ";\n\n";
 
   // exogenous variables
   os << "varexo";
-  for (int i = 0; i < (int) atoms.get_exovars().size(); i++)
-    os << " " << atoms.get_exovars()[i];
+  for (auto i : atoms.get_exovars())
+    os << " " << i;
   os << ";\n\n";
 
   // parameters
   os << "parameters";
-  for (int i = 0; i < (int) atoms.get_params().size(); i++)
-    os << " " << atoms.get_params()[i];
+  for (auto i : atoms.get_params())
+    os << " " << i;
   os << ";\n\n";
 
   // parameter values
@@ -216,9 +216,8 @@ DynareModel::check_model() const
     {
       int ft = eqs.formula(i);
       const unordered_set<int> &nuls = eqs.nulary_of_term(ft);
-      for (unordered_set<int>::const_iterator it = nuls.begin();
-           it != nuls.end(); ++it)
-        if (!atoms.is_constant(*it) && !atoms.is_named_atom(*it))
+      for (int nul : nuls)
+        if (!atoms.is_constant(nul) && !atoms.is_named_atom(nul))
           throw DynareException(__FILE__, __LINE__,
                                 "Dangling nulary term found, internal error.");
     }
@@ -268,10 +267,8 @@ DynareModel::variable_shift_map(const unordered_set<int> &a_set, int tshift,
                                 map<int, int> &s_map)
 {
   s_map.clear();
-  for (unordered_set<int>::const_iterator it = a_set.begin();
-       it != a_set.end(); ++it)
+  for (int t : a_set)
     {
-      int t = *it;
       // make shift map only for non-constants and non-parameters
       if (!atoms.is_constant(t))
         {
@@ -292,14 +289,13 @@ DynareModel::termspan(int t, int &mlead, int &mlag) const
   mlead = INT_MIN;
   mlag = INT_MAX;
   const unordered_set<int> &nul_terms = eqs.nulary_of_term(t);
-  for (unordered_set<int>::const_iterator ni = nul_terms.begin();
-       ni != nul_terms.end(); ++ni)
+  for (int nul_term : nul_terms)
     {
-      if (!atoms.is_constant(*ni)
-          && (atoms.is_type(atoms.name(*ni), DynareDynamicAtoms::endovar)
-              || atoms.is_type(atoms.name(*ni), DynareDynamicAtoms::exovar)))
+      if (!atoms.is_constant(nul_term)
+          && (atoms.is_type(atoms.name(nul_term), DynareDynamicAtoms::endovar)
+              || atoms.is_type(atoms.name(nul_term), DynareDynamicAtoms::exovar)))
         {
-          int ll = atoms.lead(*ni);
+          int ll = atoms.lead(nul_term);
           if (ll < mlag)
             mlag = ll;
           if (ll > mlead)
@@ -312,10 +308,9 @@ bool
 DynareModel::is_constant_term(int t) const
 {
   const unordered_set<int> &nul_terms = eqs.nulary_of_term(t);
-  for (unordered_set<int>::const_iterator ni = nul_terms.begin();
-       ni != nul_terms.end(); ++ni)
-    if (!atoms.is_constant(*ni)
-        && !atoms.is_type(atoms.name(*ni), DynareDynamicAtoms::param))
+  for (int nul_term : nul_terms)
+    if (!atoms.is_constant(nul_term)
+        && !atoms.is_type(atoms.name(nul_term), DynareDynamicAtoms::param))
       return false;
   return true;
 }
@@ -827,9 +822,9 @@ ModelSSWriter::write_der1(FILE *fd)
   for (int i = 0; i < model.getParser().nformulas(); i++)
     {
       const ogp::FormulaDerivatives &fder = model.getParser().derivatives(i);
-      for (unsigned int j = 0; j < eam.size(); j++)
+      for (int j : eam)
         {
-          int t = fder.derivative(ogp::FoldMultiIndex(variables.size(), 1, eam[j]));
+          int t = fder.derivative(ogp::FoldMultiIndex(variables.size(), 1, j));
           if (t > 0)
             otree.print_operation_tree(t, fd, *this);
         }
@@ -891,16 +886,14 @@ MatlabSSWriter::write_common1_preamble(FILE *fd) const
           model.getAtoms().ny(), DYNVERSION);
   // write ordering of parameters
   fprintf(fd, "\n%% params ordering\n%% =====================\n");
-  for (unsigned int ip = 0; ip < model.getAtoms().get_params().size(); ip++)
+  for (auto parname : model.getAtoms().get_params())
     {
-      const char *parname = model.getAtoms().get_params()[ip];
       fprintf(fd, "%% %s\n", parname);
     }
   // write endogenous variables
   fprintf(fd, "%%\n%% y ordering\n%% =====================\n");
-  for (unsigned int ie = 0; ie < model.getAtoms().get_endovars().size(); ie++)
+  for (auto endoname : model.getAtoms().get_endovars())
     {
-      const char *endoname = model.getAtoms().get_endovars()[ie];
       fprintf(fd, "%% %s\n", endoname);
     }
   fprintf(fd, "\n");
@@ -933,11 +926,10 @@ MatlabSSWriter::write_atom_assignment(FILE *fd) const
   // write numerical constants
   fprintf(fd, "%% numerical constants\n");
   const ogp::Constants::Tconstantmap &cmap = model.getAtoms().get_constantmap();
-  for (ogp::Constants::Tconstantmap::const_iterator it = cmap.begin();
-       it != cmap.end(); ++it)
+  for (auto it : cmap)
     {
-      format_nulary((*it).first, fd);
-      fprintf(fd, " = %12.8g;\n", (*it).second);
+      format_nulary(it.first, fd);
+      fprintf(fd, " = %12.8g;\n", it.second);
     }
   // write parameters
   fprintf(fd, "%% parameter values\n");
@@ -963,10 +955,9 @@ MatlabSSWriter::write_atom_assignment(FILE *fd) const
       try
         {
           const ogp::DynamicAtoms::Tlagmap &lmap = model.getAtoms().lagmap(exoname);
-          for (ogp::DynamicAtoms::Tlagmap::const_iterator it = lmap.begin();
-               it != lmap.end(); ++it)
+          for (auto it : lmap)
             {
-              format_nulary((*it).second, fd);
+              format_nulary(it.second, fd);
               fprintf(fd, " = 0.0; %% %s\n", exoname);
             }
         }
@@ -981,10 +972,9 @@ MatlabSSWriter::write_atom_assignment(FILE *fd) const
     {
       const char *endoname = model.getAtoms().get_endovars()[ie];
       const ogp::DynamicAtoms::Tlagmap &lmap = model.getAtoms().lagmap(endoname);
-      for (ogp::DynamicAtoms::Tlagmap::const_iterator it = lmap.begin();
-           it != lmap.end(); ++it)
+      for (auto it : lmap)
         {
-          format_nulary((*it).second, fd);
+          format_nulary(it.second, fd);
           fprintf(fd, " = y(%d); %% %s\n", ie+1, endoname);
         }
     }
@@ -1021,12 +1011,12 @@ MatlabSSWriter::write_der1_assignment(FILE *fd) const
   for (int i = 0; i < model.getParser().nformulas(); i++)
     {
       const ogp::FormulaDerivatives &fder = model.getParser().derivatives(i);
-      for (unsigned int j = 0; j < eam.size(); j++)
+      for (int j : eam)
         {
-          int tvar = variables[eam[j]];
+          int tvar = variables[j];
           const char *name = model.getAtoms().name(tvar);
           int yi = model.getAtoms().name2outer_endo(name);
-          int t = fder.derivative(ogp::FoldMultiIndex(variables.size(), 1, eam[j]));
+          int t = fder.derivative(ogp::FoldMultiIndex(variables.size(), 1, j));
           if (t != ogp::OperationTree::zero)
             {
               fprintf(fd, "out(%d,%d) = out(%d,%d) + ", i+1, yi+1, i+1, yi+1);
diff --git a/dynare++/src/forw_subst_builder.cc b/dynare++/src/forw_subst_builder.cc
index 10a97fb99..f6907c18f 100644
--- a/dynare++/src/forw_subst_builder.cc
+++ b/dynare++/src/forw_subst_builder.cc
@@ -114,20 +114,19 @@ void
 ForwSubstBuilder::unassign_gt_1_leads()
 {
   const vector<const char *> &endovars = model.atoms.get_endovars();
-  for (unsigned int i = 0; i < endovars.size(); i++)
-    unassign_gt_1_leads(endovars[i]);
+  for (auto endovar : endovars)
+    unassign_gt_1_leads(endovar);
   const vector<const char *> &exovars = model.atoms.get_exovars();
-  for (unsigned int i = 0; i < exovars.size(); i++)
-    unassign_gt_1_leads(exovars[i]);
+  for (auto exovar : exovars)
+    unassign_gt_1_leads(exovar);
 }
 
 ForwSubstBuilder::ForwSubstBuilder(const ForwSubstBuilder &b, DynareModel &m)
   : model(m)
 {
-  for (Tsubstmap::const_iterator it = b.aux_map.begin();
-       it != b.aux_map.end(); ++it)
+  for (auto it : b.aux_map)
     {
-      const char *ss = m.atoms.get_name_storage().query((*it).first);
-      aux_map.insert(Tsubstmap::value_type(ss, (*it).second));
+      const char *ss = m.atoms.get_name_storage().query(it.first);
+      aux_map.insert(Tsubstmap::value_type(ss, it.second));
     }
 }
diff --git a/dynare++/src/planner_builder.cc b/dynare++/src/planner_builder.cc
index b3d403525..fd58ffd0d 100644
--- a/dynare++/src/planner_builder.cc
+++ b/dynare++/src/planner_builder.cc
@@ -302,8 +302,8 @@ void
 PlannerBuilder::fill_yset(const ogp::NameStorage &ns,
                           const PlannerBuilder::Tvarset &yyset)
 {
-  for (Tvarset::const_iterator it = yyset.begin(); it != yyset.end(); ++it)
-    yset.insert(ns.query(*it));
+  for (auto it : yyset)
+    yset.insert(ns.query(it));
 }
 
 void
@@ -311,15 +311,13 @@ PlannerBuilder::fill_aux_map(const ogp::NameStorage &ns, const Tsubstmap &aaux_m
                              const Tsubstmap &astatic_aux_map)
 {
   // fill aux_map
-  for (Tsubstmap::const_iterator it = aaux_map.begin();
-       it != aaux_map.end(); ++it)
-    aux_map.insert(Tsubstmap::value_type(ns.query((*it).first), (*it).second));
+  for (auto it : aaux_map)
+    aux_map.insert(Tsubstmap::value_type(ns.query(it.first), it.second));
 
   // fill static_aux_map
-  for (Tsubstmap::const_iterator it = astatic_aux_map.begin();
-       it != astatic_aux_map.end(); ++it)
-    static_aux_map.insert(Tsubstmap::value_type(static_atoms.get_name_storage().query((*it).first),
-                                                (*it).second));
+  for (auto it : astatic_aux_map)
+    static_aux_map.insert(Tsubstmap::value_type(static_atoms.get_name_storage().query(it.first),
+                                                it.second));
 }
 
 MultInitSS::MultInitSS(const PlannerBuilder &pb, const Vector &pvals, Vector &yy)
@@ -380,10 +378,9 @@ MultInitSS::MultInitSS(const PlannerBuilder &pb, const Vector &pvals, Vector &yy
           if (it != old2new.end())
             {
               const ogp::AtomSubstitutions::Tshiftnameset &sset = (*it).second;
-              for (ogp::AtomSubstitutions::Tshiftnameset::const_iterator itt = sset.begin();
-                   itt != sset.end(); ++itt)
+              for (const auto & itt : sset)
                 {
-                  const char *newname = (*itt).first;
+                  const char *newname = itt.first;
                   int iouter = builder.model.atoms.name2outer_endo(newname);
                   int iy = builder.model.atoms.outer2y_endo()[iouter];
                   if (!std::isfinite(yy[iy]))
diff --git a/dynare++/sylv/cc/QuasiTriangular.cc b/dynare++/sylv/cc/QuasiTriangular.cc
index ef502b50a..e9f910f5f 100644
--- a/dynare++/sylv/cc/QuasiTriangular.cc
+++ b/dynare++/sylv/cc/QuasiTriangular.cc
@@ -100,9 +100,8 @@ Diagonal::Diagonal(double *data, const Diagonal &d)
   num_all = d.num_all;
   num_real = d.num_real;
   int d_size = d.getSize();
-  for (const_diag_iter it = d.begin(); it != d.end(); ++it)
+  for (const auto & dit : d)
     {
-      const DiagonalBlock &dit = *it;
       double *beta1 = NULL;
       double *beta2 = NULL;
       int id = dit.getIndex()*(d_size+1);
@@ -150,22 +149,22 @@ void
 Diagonal::changeBase(double *p)
 {
   int d_size = getSize();
-  for (diag_iter it = begin(); it != end(); ++it)
+  for (auto & it : *this)
     {
-      const DiagonalBlock &b = *it;
+      const DiagonalBlock &b = it;
       int jbar = b.getIndex();
       int base = d_size*jbar + jbar;
       if (b.isReal())
         {
           DiagonalBlock bnew(jbar, true, &p[base], &p[base],
                              NULL, NULL);
-          *it = bnew;
+          it = bnew;
         }
       else
         {
           DiagonalBlock bnew(jbar, false, &p[base], &p[base+d_size+1],
                              &p[base+d_size], &p[base+1]);
-          *it = bnew;
+          it = bnew;
         }
     }
 }
@@ -181,9 +180,8 @@ Diagonal::getEigenValues(Vector &eig) const
               eig.length(), 2*d_size);
       throw SYLV_MES_EXCEPTION(mes);
     }
-  for (const_diag_iter it = begin(); it != end(); ++it)
+  for (const auto & b : *this)
     {
-      const DiagonalBlock &b = *it;
       int ind = b.getIndex();
       eig[2*ind] = *(b.getAlpha());
       if (b.isReal())
@@ -307,16 +305,16 @@ void
 Diagonal::print() const
 {
   printf("Num real: %d, num complex: %d\n", getNumReal(), getNumComplex());
-  for (const_diag_iter it = begin(); it != end(); ++it)
+  for (const auto & it : *this)
     {
-      if ((*it).isReal())
+      if (it.isReal())
         {
-          printf("real: jbar=%d, alpha=%f\n", (*it).getIndex(), *((*it).getAlpha()));
+          printf("real: jbar=%d, alpha=%f\n", it.getIndex(), *(it.getAlpha()));
         }
       else
         {
           printf("complex: jbar=%d, alpha=%f, beta1=%f, beta2=%f\n",
-                 (*it).getIndex(), *((*it).getAlpha()), (*it).getBeta1(), (*it).getBeta2());
+                 it.getIndex(), *(it.getAlpha()), it.getBeta1(), it.getBeta2());
         }
     }
 }
diff --git a/dynare++/sylv/cc/Vector.cc b/dynare++/sylv/cc/Vector.cc
index d0722c5df..026553a1c 100644
--- a/dynare++/sylv/cc/Vector.cc
+++ b/dynare++/sylv/cc/Vector.cc
@@ -413,6 +413,6 @@ ConstVector::print() const
 
 ZeroPad::ZeroPad()
 {
-  for (int i = 0; i < length; i++)
-    pad[i] = 0.0;
+  for (double & i : pad)
+    i = 0.0;
 }
diff --git a/dynare++/tl/cc/equivalence.cc b/dynare++/tl/cc/equivalence.cc
index 943c2dd0c..5d9e21a2c 100644
--- a/dynare++/tl/cc/equivalence.cc
+++ b/dynare++/tl/cc/equivalence.cc
@@ -84,8 +84,8 @@ double
 OrdSequence::average() const
 {
   double res = 0;
-  for (unsigned int i = 0; i < data.size(); i++)
-    res += data[i];
+  for (int i : data)
+    res += i;
   TL_RAISE_IF(data.size() == 0,
               "Attempt to take average of empty class in OrdSequence::average");
   return res/data.size();
@@ -96,8 +96,8 @@ void
 OrdSequence::print(const char *prefix) const
 {
   printf("%s", prefix);
-  for (unsigned int i = 0; i < data.size(); i++)
-    printf("%d ", data[i]);
+  for (int i : data)
+    printf("%d ", i);
   printf("\n");
 }
 
@@ -404,8 +404,8 @@ EquivalenceBundle::EquivalenceBundle(int nmax)
 /* Destruct bundle. Just free all pointers. */
 EquivalenceBundle::~EquivalenceBundle()
 {
-  for (unsigned int i = 0; i < bundle.size(); i++)
-    delete bundle[i];
+  for (auto & i : bundle)
+    delete i;
 }
 
 /* Remember, that the first item is |EquivalenceSet(1)|. */
diff --git a/dynare++/tl/cc/fs_tensor.cc b/dynare++/tl/cc/fs_tensor.cc
index 8e7db4d62..08e389f17 100644
--- a/dynare++/tl/cc/fs_tensor.cc
+++ b/dynare++/tl/cc/fs_tensor.cc
@@ -60,11 +60,10 @@ FFSTensor::FFSTensor(const FSSparseTensor &t)
     nv(t.nvar())
 {
   zeros();
-  for (FSSparseTensor::const_iterator it = t.getMap().begin();
-       it != t.getMap().end(); ++it)
+  for (const auto & it : t.getMap())
     {
-      index ind(this, (*it).first);
-      get((*it).second.first, *ind) = (*it).second.second;
+      index ind(this, it.first);
+      get(it.second.first, *ind) = it.second.second;
     }
 }
 
diff --git a/dynare++/tl/cc/gs_tensor.cc b/dynare++/tl/cc/gs_tensor.cc
index 100788b8d..13d1ad79c 100644
--- a/dynare++/tl/cc/gs_tensor.cc
+++ b/dynare++/tl/cc/gs_tensor.cc
@@ -249,11 +249,10 @@ FGSTensor::FGSTensor(const GSSparseTensor &t)
             t.getDims().calcFoldMaxOffset(), t.dimen()), tdims(t.getDims())
 {
   zeros();
-  for (FSSparseTensor::const_iterator it = t.getMap().begin();
-       it != t.getMap().end(); ++it)
+  for (const auto & it : t.getMap())
     {
-      index ind(this, (*it).first);
-      get((*it).second.first, *ind) = (*it).second.second;
+      index ind(this, it.first);
+      get(it.second.first, *ind) = it.second.second;
     }
 }
 
diff --git a/dynare++/tl/cc/int_sequence.cc b/dynare++/tl/cc/int_sequence.cc
index 70f98e6e7..e2f3804a8 100644
--- a/dynare++/tl/cc/int_sequence.cc
+++ b/dynare++/tl/cc/int_sequence.cc
@@ -37,8 +37,8 @@ IntSequence::IntSequence(const Symmetry &sy, const vector<int> &se)
   for (int i = 0; i < length; i++)
     operator[](i) = 0;
 
-  for (unsigned int i = 0; i < se.size(); i++)
-    operator[](sy.findClass(se[i]))++;
+  for (int i : se)
+    operator[](sy.findClass(i))++;
 }
 
 /* This constructs an ordered integer sequence from the given ordered
diff --git a/dynare++/tl/cc/kron_prod.cc b/dynare++/tl/cc/kron_prod.cc
index 0b8f32d06..c286b50ca 100644
--- a/dynare++/tl/cc/kron_prod.cc
+++ b/dynare++/tl/cc/kron_prod.cc
@@ -382,8 +382,8 @@ KronProdAll::multRows(const IntSequence &irows) const
       delete row;
     }
 
-  for (unsigned int i = 0; i < to_delete.size(); i++)
-    delete to_delete[i];
+  for (auto & i : to_delete)
+    delete i;
 
   return last;
 }
diff --git a/dynare++/tl/cc/normal_moments.cc b/dynare++/tl/cc/normal_moments.cc
index 82c257a76..013d363ae 100644
--- a/dynare++/tl/cc/normal_moments.cc
+++ b/dynare++/tl/cc/normal_moments.cc
@@ -52,12 +52,11 @@ UNormalMoments::generateMoments(int maxdim, const TwoDMatrix &v)
          how the |Equivalence::apply| method works. */
       mom->zeros();
       const EquivalenceSet eset = ebundle.get(d);
-      for (EquivalenceSet::const_iterator cit = eset.begin();
-           cit != eset.end(); cit++)
+      for (const auto & cit : eset)
         {
-          if (selectEquiv(*cit))
+          if (selectEquiv(cit))
             {
-              Permutation per(*cit);
+              Permutation per(cit);
               per.inverse();
               for (Tensor::index it = kronv->begin(); it != kronv->end(); ++it)
                 {
@@ -80,10 +79,9 @@ UNormalMoments::selectEquiv(const Equivalence &e)
 {
   if (2*e.numClasses() != e.getN())
     return false;
-  for (Equivalence::const_seqit si = e.begin();
-       si != e.end(); ++si)
+  for (const auto & si : e)
     {
-      if ((*si).length() != 2)
+      if (si.length() != 2)
         return false;
     }
   return true;
@@ -94,10 +92,9 @@ UNormalMoments::selectEquiv(const Equivalence &e)
 FNormalMoments::FNormalMoments(const UNormalMoments &moms)
   : TensorContainer<FRSingleTensor>(1)
 {
-  for (UNormalMoments::const_iterator it = moms.begin();
-       it != moms.end(); ++it)
+  for (const auto & mom : moms)
     {
-      FRSingleTensor *fm = new FRSingleTensor(*((*it).second));
+      FRSingleTensor *fm = new FRSingleTensor(*(mom.second));
       insert(fm);
     }
 }
diff --git a/dynare++/tl/cc/permutation.cc b/dynare++/tl/cc/permutation.cc
index 43d72372c..d94cceedd 100644
--- a/dynare++/tl/cc/permutation.cc
+++ b/dynare++/tl/cc/permutation.cc
@@ -131,8 +131,8 @@ PermutationBundle::PermutationBundle(int nmax)
 
 PermutationBundle::~PermutationBundle()
 {
-  for (unsigned int i = 0; i < bundle.size(); i++)
-    delete bundle[i];
+  for (auto & i : bundle)
+    delete i;
 }
 
 const PermutationSet &
diff --git a/dynare++/tl/cc/ps_tensor.cc b/dynare++/tl/cc/ps_tensor.cc
index 3d191d0b3..c9a57b8a0 100644
--- a/dynare++/tl/cc/ps_tensor.cc
+++ b/dynare++/tl/cc/ps_tensor.cc
@@ -228,10 +228,10 @@ UPSTensor::fillFromSparseTwo(const FSSparseTensor &t, const IntSequence &ss,
           IntSequence c((*run).first);
           c.add(-1, lb_srt);
           unsort.apply(c);
-          for (unsigned int i = 0; i < pp.size(); i++)
+          for (auto & i : pp)
             {
               IntSequence cp(coor.size());
-              pp[i]->apply(c, cp);
+              i->apply(c, cp);
               Tensor::index ind(this, cp);
               TL_RAISE_IF(*ind < 0 || *ind >= ncols(),
                           "Internal error in slicing constructor of UPSTensor");
diff --git a/dynare++/tl/cc/pyramid_prod.cc b/dynare++/tl/cc/pyramid_prod.cc
index 651349378..a7ba4ba14 100644
--- a/dynare++/tl/cc/pyramid_prod.cc
+++ b/dynare++/tl/cc/pyramid_prod.cc
@@ -26,14 +26,13 @@ USubTensor::USubTensor(const TensorDimens &bdims,
               "Tensor has not full symmetry in USubTensor()");
   const EquivalenceSet &eset = cont.getEqBundle().get(bdims.dimen());
   zeros();
-  for (EquivalenceSet::const_iterator it = eset.begin();
-       it != eset.end(); ++it)
+  for (const auto & it : eset)
     {
-      if ((*it).numClasses() == hdims.dimen())
+      if (it.numClasses() == hdims.dimen())
         {
-          Permutation per(*it);
+          Permutation per(it);
           vector<const FGSTensor *> ts
-            = cont.fetchTensors(bdims.getSym(), *it);
+            = cont.fetchTensors(bdims.getSym(), it);
           for (int i = 0; i < (int) lst.size(); i++)
             {
               IntSequence perindex(lst[i].size());
@@ -64,12 +63,12 @@ USubTensor::addKronColumn(int i, const vector<const FGSTensor *> &ts,
 {
   vector<ConstVector> tmpcols;
   int lastdim = 0;
-  for (unsigned int j = 0; j < ts.size(); j++)
+  for (auto t : ts)
     {
-      IntSequence ind(pindex, lastdim, lastdim+ts[j]->dimen());
-      lastdim += ts[j]->dimen();
-      index in(ts[j], ind);
-      tmpcols.push_back(ConstVector(*(ts[j]), *in));
+      IntSequence ind(pindex, lastdim, lastdim+t->dimen());
+      lastdim += t->dimen();
+      index in(t, ind);
+      tmpcols.push_back(ConstVector(*t, *in));
     }
 
   URSingleTensor kronmult(tmpcols);
diff --git a/dynare++/tl/cc/stack_container.cc b/dynare++/tl/cc/stack_container.cc
index 126069df1..ee3f192dd 100644
--- a/dynare++/tl/cc/stack_container.cc
+++ b/dynare++/tl/cc/stack_container.cc
@@ -111,12 +111,11 @@ WorkerFoldMAASparse1::operator()()
       const Permutation &per = pset.get(iper);
       IntSequence percoor(coor.size());
       per.apply(coor, percoor);
-      for (EquivalenceSet::const_iterator it = eset.begin();
-           it != eset.end(); ++it)
+      for (const auto & it : eset)
         {
-          if ((*it).numClasses() == t.dimen())
+          if (it.numClasses() == t.dimen())
             {
-              StackProduct<FGSTensor> sp(cont, *it, out.getSym());
+              StackProduct<FGSTensor> sp(cont, it, out.getSym());
               if (!sp.isZero(percoor))
                 {
                   KronProdStack<FGSTensor> kp(sp, percoor);
@@ -124,7 +123,7 @@ WorkerFoldMAASparse1::operator()()
                   const Permutation &oper = kp.getPer();
                   if (Permutation(oper, per) == iden)
                     {
-                      FPSTensor fps(out.getDims(), *it, slice, kp);
+                      FPSTensor fps(out.getDims(), it, slice, kp);
                       {
                         SYNCHRO syn(&out, "WorkerUnfoldMAASparse1");
                         fps.addTo(out);
@@ -226,12 +225,11 @@ FoldedStackContainer::multAndAddSparse3(const FSSparseTensor &t,
       Vector outcol(out, *run);
       FRSingleTensor sumcol(t.nvar(), t.dimen());
       sumcol.zeros();
-      for (EquivalenceSet::const_iterator it = eset.begin();
-           it != eset.end(); ++it)
+      for (const auto & it : eset)
         {
-          if ((*it).numClasses() == t.dimen())
+          if (it.numClasses() == t.dimen())
             {
-              StackProduct<FGSTensor> sp(*this, *it, out.getSym());
+              StackProduct<FGSTensor> sp(*this, it, out.getSym());
               IrregTensorHeader header(sp, run.getCoor());
               IrregTensor irten(header);
               irten.addTo(sumcol);
@@ -308,18 +306,17 @@ FoldedStackContainer::multAndAddStacks(const IntSequence &coor,
         {
           Permutation sort_per(ui.getCoor());
           sort_per.inverse();
-          for (EquivalenceSet::const_iterator it = eset.begin();
-               it != eset.end(); ++it)
+          for (const auto & it : eset)
             {
-              if ((*it).numClasses() == g.dimen())
+              if (it.numClasses() == g.dimen())
                 {
-                  StackProduct<FGSTensor> sp(*this, *it, sort_per, out.getSym());
+                  StackProduct<FGSTensor> sp(*this, it, sort_per, out.getSym());
                   if (!sp.isZero(coor))
                     {
                       KronProdStack<FGSTensor> kp(sp, coor);
                       if (ug.getSym().isFull())
                         kp.optimizeOrder();
-                      FPSTensor fps(out.getDims(), *it, sort_per, ug, kp);
+                      FPSTensor fps(out.getDims(), it, sort_per, ug, kp);
                       {
                         SYNCHRO syn(ad, "multAndAddStacks");
                         fps.addTo(out);
@@ -352,16 +349,15 @@ FoldedStackContainer::multAndAddStacks(const IntSequence &coor,
         {
           Permutation sort_per(ui.getCoor());
           sort_per.inverse();
-          for (EquivalenceSet::const_iterator it = eset.begin();
-               it != eset.end(); ++it)
+          for (const auto & it : eset)
             {
-              if ((*it).numClasses() == g.dimen())
+              if (it.numClasses() == g.dimen())
                 {
-                  StackProduct<FGSTensor> sp(*this, *it, sort_per, out.getSym());
+                  StackProduct<FGSTensor> sp(*this, it, sort_per, out.getSym());
                   if (!sp.isZero(coor))
                     {
                       KronProdStack<FGSTensor> kp(sp, coor);
-                      FPSTensor fps(out.getDims(), *it, sort_per, g, kp);
+                      FPSTensor fps(out.getDims(), it, sort_per, g, kp);
                       {
                         SYNCHRO syn(ad, "multAndAddStacks");
                         fps.addTo(out);
@@ -507,12 +503,11 @@ WorkerUnfoldMAASparse1::operator()()
       const Permutation &per = pset.get(iper);
       IntSequence percoor(coor.size());
       per.apply(coor, percoor);
-      for (EquivalenceSet::const_iterator it = eset.begin();
-           it != eset.end(); ++it)
+      for (const auto & it : eset)
         {
-          if ((*it).numClasses() == t.dimen())
+          if (it.numClasses() == t.dimen())
             {
-              StackProduct<UGSTensor> sp(cont, *it, out.getSym());
+              StackProduct<UGSTensor> sp(cont, it, out.getSym());
               if (!sp.isZero(percoor))
                 {
                   KronProdStack<UGSTensor> kp(sp, percoor);
@@ -520,7 +515,7 @@ WorkerUnfoldMAASparse1::operator()()
                   const Permutation &oper = kp.getPer();
                   if (Permutation(oper, per) == iden)
                     {
-                      UPSTensor ups(out.getDims(), *it, slice, kp);
+                      UPSTensor ups(out.getDims(), it, slice, kp);
                       {
                         SYNCHRO syn(&out, "WorkerUnfoldMAASparse1");
                         ups.addTo(out);
@@ -638,18 +633,17 @@ UnfoldedStackContainer::multAndAddStacks(const IntSequence &fi,
         {
           Permutation sort_per(ui.getCoor());
           sort_per.inverse();
-          for (EquivalenceSet::const_iterator it = eset.begin();
-               it != eset.end(); ++it)
+          for (const auto & it : eset)
             {
-              if ((*it).numClasses() == g.dimen())
+              if (it.numClasses() == g.dimen())
                 {
-                  StackProduct<UGSTensor> sp(*this, *it, sort_per, out.getSym());
+                  StackProduct<UGSTensor> sp(*this, it, sort_per, out.getSym());
                   if (!sp.isZero(fi))
                     {
                       KronProdStack<UGSTensor> kp(sp, fi);
                       if (g.getSym().isFull())
                         kp.optimizeOrder();
-                      UPSTensor ups(out.getDims(), *it, sort_per, g, kp);
+                      UPSTensor ups(out.getDims(), it, sort_per, g, kp);
                       {
                         SYNCHRO syn(ad, "multAndAddStacks");
                         ups.addTo(out);
diff --git a/dynare++/tl/cc/symmetry.cc b/dynare++/tl/cc/symmetry.cc
index 8aafd2161..e30577f6c 100644
--- a/dynare++/tl/cc/symmetry.cc
+++ b/dynare++/tl/cc/symmetry.cc
@@ -116,9 +116,9 @@ symiterator::operator++()
 
 InducedSymmetries::InducedSymmetries(const Equivalence &e, const Symmetry &s)
 {
-  for (Equivalence::const_seqit i = e.begin(); i != e.end(); ++i)
+  for (const auto & i : e)
     {
-      push_back(Symmetry(s, *i));
+      push_back(Symmetry(s, i));
     }
 }
 
diff --git a/dynare++/tl/cc/t_container.cc b/dynare++/tl/cc/t_container.cc
index 007f958e8..55aed4623 100644
--- a/dynare++/tl/cc/t_container.cc
+++ b/dynare++/tl/cc/t_container.cc
@@ -11,10 +11,9 @@ const int FGSContainer::num_one_time = 10;
 UGSContainer::UGSContainer(const FGSContainer &c)
   : TensorContainer<UGSTensor>(c.num())
 {
-  for (FGSContainer::const_iterator it = c.begin();
-       it != c.end(); ++it)
+  for (const auto & it : c)
     {
-      UGSTensor *unfolded = new UGSTensor(*((*it).second));
+      UGSTensor *unfolded = new UGSTensor(*(it.second));
       insert(unfolded);
     }
 }
@@ -38,18 +37,17 @@ UGSContainer::multAndAdd(const UGSTensor &t, UGSTensor &out) const
   int k = out.dimen();
   const EquivalenceSet &eset = ebundle.get(k);
 
-  for (EquivalenceSet::const_iterator it = eset.begin();
-       it != eset.end(); ++it)
+  for (const auto & it : eset)
     {
-      if ((*it).numClasses() == l)
+      if (it.numClasses() == l)
         {
           vector<const UGSTensor *> ts
-            = fetchTensors(out.getSym(), *it);
+            = fetchTensors(out.getSym(), it);
           KronProdAllOptim kp(l);
           for (int i = 0; i < l; i++)
             kp.setMat(i, *(ts[i]));
           kp.optimizeOrder();
-          UPSTensor ups(out.getDims(), *it, t, kp);
+          UPSTensor ups(out.getDims(), it, t, kp);
           ups.addTo(out);
         }
     }
@@ -59,10 +57,9 @@ UGSContainer::multAndAdd(const UGSTensor &t, UGSTensor &out) const
 FGSContainer::FGSContainer(const UGSContainer &c)
   : TensorContainer<FGSTensor>(c.num())
 {
-  for (UGSContainer::const_iterator it = c.begin();
-       it != c.end(); ++it)
+  for (const auto & it : c)
     {
-      FGSTensor *folded = new FGSTensor(*((*it).second));
+      FGSTensor *folded = new FGSTensor(*(it.second));
       insert(folded);
     }
 }
@@ -88,18 +85,17 @@ FGSContainer::multAndAdd(const UGSTensor &t, FGSTensor &out) const
   int k = out.dimen();
   const EquivalenceSet &eset = ebundle.get(k);
 
-  for (EquivalenceSet::const_iterator it = eset.begin();
-       it != eset.end(); ++it)
+  for (const auto & it : eset)
     {
-      if ((*it).numClasses() == l)
+      if (it.numClasses() == l)
         {
           vector<const FGSTensor *> ts
-            = fetchTensors(out.getSym(), *it);
+            = fetchTensors(out.getSym(), it);
           KronProdAllOptim kp(l);
           for (int i = 0; i < l; i++)
             kp.setMat(i, *(ts[i]));
           kp.optimizeOrder();
-          FPSTensor fps(out.getDims(), *it, t, kp);
+          FPSTensor fps(out.getDims(), it, t, kp);
           fps.addTo(out);
         }
     }
diff --git a/dynare++/tl/cc/t_polynomial.cc b/dynare++/tl/cc/t_polynomial.cc
index 2ab0e7534..f2ff1fbf7 100644
--- a/dynare++/tl/cc/t_polynomial.cc
+++ b/dynare++/tl/cc/t_polynomial.cc
@@ -50,19 +50,17 @@ PowerProvider::~PowerProvider()
 UTensorPolynomial::UTensorPolynomial(const FTensorPolynomial &fp)
   : TensorPolynomial<UFSTensor, UGSTensor, URSingleTensor>(fp.nrows(), fp.nvars())
 {
-  for (FTensorPolynomial::const_iterator it = fp.begin();
-       it != fp.end(); ++it)
+  for (const auto & it : fp)
     {
-      insert(new UFSTensor(*((*it).second)));
+      insert(new UFSTensor(*(it.second)));
     }
 }
 
 FTensorPolynomial::FTensorPolynomial(const UTensorPolynomial &up)
   : TensorPolynomial<FFSTensor, FGSTensor, FRSingleTensor>(up.nrows(), up.nvars())
 {
-  for (UTensorPolynomial::const_iterator it = up.begin();
-       it != up.end(); ++it)
+  for (const auto & it : up)
     {
-      insert(new FFSTensor(*((*it).second)));
+      insert(new FFSTensor(*(it.second)));
     }
 }
diff --git a/dynare++/utils/cc/pascal_triangle.cc b/dynare++/utils/cc/pascal_triangle.cc
index c64db047d..21e1d36dd 100644
--- a/dynare++/utils/cc/pascal_triangle.cc
+++ b/dynare++/utils/cc/pascal_triangle.cc
@@ -96,6 +96,6 @@ PascalTriangle::noverk(int n, int k)
 void
 PascalTriangle::print() const
 {
-  for (unsigned int i = 0; i < tr.size(); i++)
-    tr[i].print();
+  for (const auto & i : tr)
+    i.print();
 }