diff --git a/dynare++/integ/cc/product.cc b/dynare++/integ/cc/product.cc
index d1f368c1a..aaf35362e 100644
--- a/dynare++/integ/cc/product.cc
+++ b/dynare++/integ/cc/product.cc
@@ -10,7 +10,7 @@
 
 prodpit::prodpit(const ProductQuadrature &q, int j0, int l)
   : prodq(q), level(l), npoints(q.uquad.numPoints(l)),
-    jseq{q.dimen(), 0},
+    jseq(q.dimen(), 0),
     end_flag(false),
     sig{q.dimen()},
     p{q.dimen()}
diff --git a/dynare++/integ/cc/smolyak.cc b/dynare++/integ/cc/smolyak.cc
index c155e64e7..6fec35e9c 100644
--- a/dynare++/integ/cc/smolyak.cc
+++ b/dynare++/integ/cc/smolyak.cc
@@ -12,7 +12,7 @@
 
 smolpit::smolpit(const SmolyakQuadrature &q, unsigned int isum)
   : smolq(q), isummand(isum),
-    jseq{q.dimen(), 0},
+    jseq(q.dimen(), 0),
     sig{q.dimen()},
     p{q.dimen()}
 {
diff --git a/dynare++/integ/testing/tests.cc b/dynare++/integ/testing/tests.cc
index a91e882b9..3d2ecc8a3 100644
--- a/dynare++/integ/testing/tests.cc
+++ b/dynare++/integ/testing/tests.cc
@@ -256,7 +256,7 @@ TestRunnable::smolyak_normal_moments(const GeneralMatrix &m, int imom, int level
 
   // check against theoretical moments
   UNormalMoments moments(imom, msq);
-  smol_out.add(-1.0, (moments.get(Symmetry(imom)))->getData());
+  smol_out.add(-1.0, (moments.get(Symmetry{imom}))->getData());
   std::cout << "\tError:                         " << std::setw(16) << std::setprecision(12) << smol_out.getMax() << std::endl;
   return smol_out.getMax() < 1.e-7;
 }
@@ -285,7 +285,7 @@ TestRunnable::product_normal_moments(const GeneralMatrix &m, int imom, int level
 
   // check against theoretical moments
   UNormalMoments moments(imom, msq);
-  prod_out.add(-1.0, (moments.get(Symmetry(imom)))->getData());
+  prod_out.add(-1.0, (moments.get(Symmetry{imom}))->getData());
   std::cout << "\tError:                         " << std::setw(16) << std::setprecision(12) << prod_out.getMax() << std::endl;
   return prod_out.getMax() < 1.e-7;
 }
diff --git a/dynare++/kord/approximation.cc b/dynare++/kord/approximation.cc
index ddf9a5058..3ed996d4b 100644
--- a/dynare++/kord/approximation.cc
+++ b/dynare++/kord/approximation.cc
@@ -80,7 +80,7 @@ Approximation::approxAtSteady()
 {
   model.calcDerivativesAtSteady();
   FirstOrder fo(model.nstat(), model.npred(), model.nboth(), model.nforw(),
-                model.nexog(), *(model.getModelDerivatives().get(Symmetry(1))),
+                model.nexog(), *(model.getModelDerivatives().get(Symmetry{1})),
                 journal, qz_criterium);
   KORD_RAISE_IF_X(!fo.isStable(),
                   "The model is not Blanchard-Kahn stable",
@@ -293,21 +293,21 @@ Approximation::calcStochShift(Vector &out, double at_sigma) const
     {
       if (KOrder::is_even(d))
         {
-          Symmetry sym(0, d, 0, 0);
+          Symmetry sym{0, d, 0, 0};
 
           // calculate $F_{u'^d}$ via |ZAuxContainer|
           FGSTensor *ten = new FGSTensor(ypart.ny(), TensorDimens(sym, nvs));
           ten->zeros();
           for (int l = 1; l <= d; l++)
             {
-              const FSSparseTensor *f = model.getModelDerivatives().get(Symmetry(l));
+              const FSSparseTensor *f = model.getModelDerivatives().get(Symmetry{l});
               zaux.multAndAdd(*f, *ten);
             }
 
           // multiply with shocks and add to result
-          FGSTensor *tmp = new FGSTensor(ypart.ny(), TensorDimens(Symmetry(0, 0, 0, 0), nvs));
+          FGSTensor *tmp = new FGSTensor(ypart.ny(), TensorDimens(Symmetry{0, 0, 0, 0}, nvs));
           tmp->zeros();
-          ten->contractAndAdd(1, *tmp, *(mom.get(Symmetry(d))));
+          ten->contractAndAdd(1, *tmp, *(mom.get(Symmetry{d})));
 
           out.add(pow(at_sigma, d)/dfac, tmp->getData());
           delete ten;
@@ -360,8 +360,8 @@ Approximation::check(double at_sigma) const
 TwoDMatrix *
 Approximation::calcYCov() const
 {
-  const TwoDMatrix &gy = *(rule_ders->get(Symmetry(1, 0, 0, 0)));
-  const TwoDMatrix &gu = *(rule_ders->get(Symmetry(0, 1, 0, 0)));
+  const TwoDMatrix &gy = *(rule_ders->get(Symmetry{1, 0, 0, 0}));
+  const TwoDMatrix &gu = *(rule_ders->get(Symmetry{0, 1, 0, 0}));
   TwoDMatrix G(model.numeq(), model.numeq());
   G.zeros();
   G.place(gy, 0, model.nstat());
diff --git a/dynare++/kord/decision_rule.hh b/dynare++/kord/decision_rule.hh
index 359447fb8..749f7e175 100644
--- a/dynare++/kord/decision_rule.hh
+++ b/dynare++/kord/decision_rule.hh
@@ -217,11 +217,11 @@ DecisionRuleImpl<t>::fillTensors(const _Tg &g, double sigma)
           int j = d-i;
           int kfact = 1;
           _Ttensor tmp(ypart.ny(),
-                       TensorDimens(Symmetry(i, j), tns));
+                       TensorDimens(Symmetry{i, j}, tns));
           tmp.zeros();
           for (int k = 0; k+d <= g.getMaxDim(); k++, kfact *= k)
             {
-              Symmetry sym(i, j, 0, k);
+              Symmetry sym{i, j, 0, k};
               if (g.check(sym))
                 {
                   double mult = pow(sigma, k)/dfact/kfact;
@@ -539,7 +539,7 @@ DRFixPoint<t>::DRFixPoint(const _Tg &g, const PartitionY &yp,
   fillTensors(g, sigma);
   _Tparent yspol(ypart.nstat, ypart.nys(), *this);
   bigf = new _Tparent((const _Tparent &) yspol);
-  _Ttensym *frst = bigf->get(Symmetry(1));
+  _Ttensym *frst = bigf->get(Symmetry{1});
   for (int i = 0; i < ypart.nys(); i++)
     frst->get(i, i) = frst->get(i, i) - 1;
   bigfder = new _Tparent(*bigf, 0);
@@ -572,9 +572,9 @@ DRFixPoint<t>::fillTensors(const _Tg &g, double sigma)
       int kfact = 1;
       for (int k = 0; d+k <= g.getMaxDim(); k++, kfact *= k)
         {
-          if (g.check(Symmetry(d, 0, 0, k)))
+          if (g.check(Symmetry{d, 0, 0, k}))
             {
-              const _Ttensor *ten = g.get(Symmetry(d, 0, 0, k));
+              const _Ttensor *ten = g.get(Symmetry{d, 0, 0, k});
               double mult = pow(sigma, k)/dfact/kfact;
               g_yd->add(mult, *ten);
             }
diff --git a/dynare++/kord/first_order.hh b/dynare++/kord/first_order.hh
index 60b53d446..c396570dc 100644
--- a/dynare++/kord/first_order.hh
+++ b/dynare++/kord/first_order.hh
@@ -72,10 +72,10 @@ public:
   {
     IntSequence nvs(4);
     nvs[0] = fo.ypart.nys(); nvs[1] = fo.nu; nvs[2] = fo.nu; nvs[3] = 1;
-    _Ttensor *ten = new _Ttensor(fo.ypart.ny(), TensorDimens(Symmetry(1, 0, 0, 0), nvs));
+    _Ttensor *ten = new _Ttensor(fo.ypart.ny(), TensorDimens(Symmetry{1, 0, 0, 0}, nvs));
     ten->zeros(); ten->add(1.0, fo.gy);
     this->insert(ten);
-    ten = new _Ttensor(fo.ypart.ny(), TensorDimens(Symmetry(0, 1, 0, 0), nvs));
+    ten = new _Ttensor(fo.ypart.ny(), TensorDimens(Symmetry{0, 1, 0, 0}, nvs));
     ten->zeros(); ten->add(1.0, fo.gu);
     this->insert(ten);
   }
diff --git a/dynare++/kord/global_check.cc b/dynare++/kord/global_check.cc
index 94b957e82..08fdca3bc 100644
--- a/dynare++/kord/global_check.cc
+++ b/dynare++/kord/global_check.cc
@@ -95,9 +95,9 @@ ResidFunction::setYU(const ConstVector &ys, const ConstVector &xx)
   hss = new FTensorPolynomial(dr_ss, ytmp_star);
   ConstVector ysteady_ss(dr.c->getSteady(), model->nstat()+model->npred(),
                          model->nboth()+model->nforw());
-  if (hss->check(Symmetry(0)))
+  if (hss->check(Symmetry{0}))
     {
-      hss->get(Symmetry(0))->getData().add(1.0, ysteady_ss);
+      hss->get(Symmetry{0})->getData().add(1.0, ysteady_ss);
     }
   else
     {
diff --git a/dynare++/kord/korder.cc b/dynare++/kord/korder.cc
index 89893f7bc..7f4be316a 100644
--- a/dynare++/kord/korder.cc
+++ b/dynare++/kord/korder.cc
@@ -239,9 +239,9 @@ KOrder::KOrder(int num_stat, int num_pred, int num_both, int num_forw,
     _uGstack(&_ugs, ypart.nys(), nu),
     _fGstack(&_fgs, ypart.nys(), nu),
     _um(maxk, v), _fm(_um), f(fcont),
-    matA(*(f.get(Symmetry(1))), _uZstack.getStackSizes(), gy, ypart),
-    matS(*(f.get(Symmetry(1))), _uZstack.getStackSizes(), gy, ypart),
-    matB(*(f.get(Symmetry(1))), _uZstack.getStackSizes()),
+    matA(*(f.get(Symmetry{1})), _uZstack.getStackSizes(), gy, ypart),
+    matS(*(f.get(Symmetry{1})), _uZstack.getStackSizes(), gy, ypart),
+    matB(*(f.get(Symmetry{1})), _uZstack.getStackSizes()),
     journal(jr)
 {
   KORD_RAISE_IF(gy.ncols() != ypart.nys(),
@@ -265,20 +265,20 @@ KOrder::KOrder(int num_stat, int num_pred, int num_both, int num_forw,
   // put $g_y$ and $g_u$ to the container
   /* Note that $g_\sigma$ is zero by the nature and we do not insert it to
      the container. We insert a new physical copies. */
-  UGSTensor *tgy = new UGSTensor(ny, TensorDimens(Symmetry(1, 0, 0, 0), nvs));
+  UGSTensor *tgy = new UGSTensor(ny, TensorDimens(Symmetry{1, 0, 0, 0}, nvs));
   tgy->getData() = gy.getData();
   insertDerivative<unfold>(tgy);
-  UGSTensor *tgu = new UGSTensor(ny, TensorDimens(Symmetry(0, 1, 0, 0), nvs));
+  UGSTensor *tgu = new UGSTensor(ny, TensorDimens(Symmetry{0, 1, 0, 0}, nvs));
   tgu->getData() = gu.getData();
   insertDerivative<unfold>(tgu);
 
   // put $G_y$, $G_u$ and $G_{u'}$ to the container
   /* Also note that since $g_\sigma$ is zero, so $G_\sigma$. */
-  UGSTensor *tGy = faaDiBrunoG<unfold>(Symmetry(1, 0, 0, 0));
+  UGSTensor *tGy = faaDiBrunoG<unfold>(Symmetry{1, 0, 0, 0});
   G<unfold>().insert(tGy);
-  UGSTensor *tGu = faaDiBrunoG<unfold>(Symmetry(0, 1, 0, 0));
+  UGSTensor *tGu = faaDiBrunoG<unfold>(Symmetry{0, 1, 0, 0});
   G<unfold>().insert(tGu);
-  UGSTensor *tGup = faaDiBrunoG<unfold>(Symmetry(0, 0, 1, 0));
+  UGSTensor *tGup = faaDiBrunoG<unfold>(Symmetry{0, 0, 1, 0});
   G<unfold>().insert(tGup);
 }
 
@@ -306,7 +306,7 @@ KOrder::sylvesterSolve<KOrder::unfold>(ctraits<unfold>::Ttensor &der) const
     {
       KORD_RAISE_IF(!der.isFinite(),
                     "RHS of Sylverster is not finite");
-      TwoDMatrix gs_y(*(gs<unfold>().get(Symmetry(1, 0, 0, 0))));
+      TwoDMatrix gs_y(*(gs<unfold>().get(Symmetry{1, 0, 0, 0})));
       GeneralSylvester sylv(der.getSym()[0], ny, ypart.nys(),
                             ypart.nstat+ypart.npred,
                             matA.getData(), matB.getData(),
diff --git a/dynare++/kord/korder.hh b/dynare++/kord/korder.hh
index ab9201554..7c504c00d 100644
--- a/dynare++/kord/korder.hh
+++ b/dynare++/kord/korder.hh
@@ -461,7 +461,7 @@ template<int t>
 void
 KOrder::recover_y(int i)
 {
-  Symmetry sym(i, 0, 0, 0);
+  Symmetry sym{i, 0, 0, 0};
   JournalRecordPair pa(journal);
   pa << "Recovering symmetry " << sym << endrec;
 
@@ -475,7 +475,7 @@ KOrder::recover_y(int i)
 
   insertDerivative<t>(g_yi);
 
-  _Ttensor *gss_y = gss<t>().get(Symmetry(1, 0, 0, 0));
+  _Ttensor *gss_y = gss<t>().get(Symmetry{1, 0, 0, 0});
   gs<t>().multAndAdd(*gss_y, *G_yi);
   _Ttensor *gss_yi = gss<t>().get(sym);
   gs<t>().multAndAdd(*gss_yi, *G_yi);
@@ -496,7 +496,7 @@ template <int t>
 void
 KOrder::recover_yu(int i, int j)
 {
-  Symmetry sym(i, j, 0, 0);
+  Symmetry sym{i, j, 0, 0};
   JournalRecordPair pa(journal);
   pa << "Recovering symmetry " << sym << endrec;
 
@@ -508,7 +508,7 @@ KOrder::recover_yu(int i, int j)
   matA.multInv(*g_yiuj);
   insertDerivative<t>(g_yiuj);
 
-  gs<t>().multAndAdd(*(gss<t>().get(Symmetry(1, 0, 0, 0))), *G_yiuj);
+  gs<t>().multAndAdd(*(gss<t>().get(Symmetry{1, 0, 0, 0})), *G_yiuj);
 }
 
 /* Here we solve
@@ -535,7 +535,7 @@ template <int t>
 void
 KOrder::recover_ys(int i, int j)
 {
-  Symmetry sym(i, 0, 0, j);
+  Symmetry sym{i, 0, 0, j};
   JournalRecordPair pa(journal);
   pa << "Recovering symmetry " << sym << endrec;
 
@@ -595,7 +595,7 @@ template <int t>
 void
 KOrder::recover_yus(int i, int j, int k)
 {
-  Symmetry sym(i, j, 0, k);
+  Symmetry sym{i, j, 0, k};
   JournalRecordPair pa(journal);
   pa << "Recovering symmetry " << sym << endrec;
 
@@ -662,7 +662,7 @@ template <int t>
 void
 KOrder::recover_s(int i)
 {
-  Symmetry sym(0, 0, 0, i);
+  Symmetry sym{0, 0, 0, i};
   JournalRecordPair pa(journal);
   pa << "Recovering symmetry " << sym << endrec;
 
@@ -710,7 +710,7 @@ KOrder::fillG(int i, int j, int k)
     {
       if (is_even(k-m))
         {
-          _Ttensor *G_yiujupms = faaDiBrunoG<t>(Symmetry(i, j, m, k-m));
+          _Ttensor *G_yiujupms = faaDiBrunoG<t>(Symmetry{i, j, m, k-m});
           G<t>().insert(G_yiujupms);
         }
     }
@@ -727,12 +727,12 @@ template <int t>
 _Ttensor *
 KOrder::calcD_ijk(int i, int j, int k) const
 {
-  _Ttensor *res =       new _Ttensor(ny, TensorDimens(Symmetry(i, j, 0, 0), nvs));
+  _Ttensor *res =       new _Ttensor(ny, TensorDimens(Symmetry{i, j, 0, 0}, nvs));
   res->zeros();
   if (is_even(k))
     {
-      _Ttensor *tmp = faaDiBrunoZ<t>(Symmetry(i, j, k, 0));
-      tmp->contractAndAdd(2, *res, *(m<t>().get(Symmetry(k))));
+      _Ttensor *tmp = faaDiBrunoZ<t>(Symmetry{i, j, k, 0});
+      tmp->contractAndAdd(2, *res, *(m<t>().get(Symmetry{k})));
       delete tmp;
     }
   return res;
@@ -749,13 +749,13 @@ template <int t>
 _Ttensor *
 KOrder::calcE_ijk(int i, int j, int k) const
 {
-  _Ttensor *res = new _Ttensor(ny, TensorDimens(Symmetry(i, j, 0, 0), nvs));
+  _Ttensor *res = new _Ttensor(ny, TensorDimens(Symmetry{i, j, 0, 0}, nvs));
   res->zeros();
   for (int n = 2; n <= k-1; n += 2)
     {
-      _Ttensor *tmp = faaDiBrunoZ<t>(Symmetry(i, j, n, k-n));
+      _Ttensor *tmp = faaDiBrunoZ<t>(Symmetry{i, j, n, k-n});
       tmp->mult((double) (Tensor::noverk(k, n)));
-      tmp->contractAndAdd(2, *res, *(m<t>().get(Symmetry(n))));
+      tmp->contractAndAdd(2, *res, *(m<t>().get(Symmetry{n})));
       delete tmp;
     }
   return res;
@@ -861,7 +861,7 @@ KOrder::check(int dim) const
   // check for $F_{y^iu^j}=0
   for (int i = 0; i <= dim; i++)
     {
-      Symmetry sym(dim-i, i, 0,  0);
+      Symmetry sym{dim-i, i, 0,  0};
       _Ttensor *r = faaDiBrunoZ<t>(sym);
       double err = r->getData().getMax();
       JournalRecord(journal) << "\terror for symmetry " << sym << "\tis " << err << endrec;
@@ -879,7 +879,7 @@ KOrder::check(int dim) const
       int k = (*si)[2];
       if (i+j > 0 && k > 0)
         {
-          Symmetry sym(i, j, 0, k);
+          Symmetry sym{i, j, 0, k};
           _Ttensor *r = faaDiBrunoZ<t>(sym);
           _Ttensor *D_ijk = calcD_ijk<t>(i, j, k);
           r->add(1.0, *D_ijk);
@@ -894,7 +894,7 @@ KOrder::check(int dim) const
     }
 
   // check for $F_{\sigma^i}+D_i+E_i=0
-  _Ttensor *r = faaDiBrunoZ<t>(Symmetry(0, 0, 0, dim));
+  _Ttensor *r = faaDiBrunoZ<t>(Symmetry{0, 0, 0, dim});
   _Ttensor *D_k = calcD_k<t>(dim);
   r->add(1.0, *D_k);
   delete D_k;
@@ -902,7 +902,7 @@ KOrder::check(int dim) const
   r->add(1.0, *E_k);
   delete E_k;
   double err = r->getData().getMax();
-  Symmetry sym(0, 0, 0, dim);
+  Symmetry sym{0, 0, 0, dim};
   JournalRecord(journal) << "\terror for symmetry " << sym << "\tis " << err << endrec;
   if (err > maxerror)
     maxerror = err;
diff --git a/dynare++/kord/korder_stoch.cc b/dynare++/kord/korder_stoch.cc
index ae0880c76..bf9b68fd9 100644
--- a/dynare++/kord/korder_stoch.cc
+++ b/dynare++/kord/korder_stoch.cc
@@ -35,7 +35,7 @@ KOrderStoch::KOrderStoch(const PartitionY &yp, int nu,
     _uGstack(&_ugs, ypart.nys(), nu),
     _fGstack(&_fgs, ypart.nys(), nu),
     f(fcont),
-    matA(*(fcont.get(Symmetry(1))), _uZstack.getStackSizes(), *(hh.get(Symmetry(1, 0, 0, 0))),
+    matA(*(fcont.get(Symmetry{1})), _uZstack.getStackSizes(), *(hh.get(Symmetry{1, 0, 0, 0})),
          ypart)
 {
   nvs[0] = ypart.nys();
@@ -56,7 +56,7 @@ KOrderStoch::KOrderStoch(const PartitionY &yp, int nu,
     _uGstack(&_ugs, ypart.nys(), nu),
     _fGstack(&_fgs, ypart.nys(), nu),
     f(fcont),
-    matA(*(fcont.get(Symmetry(1))), _uZstack.getStackSizes(), *(hh.get(Symmetry(1, 0, 0, 0))),
+    matA(*(fcont.get(Symmetry{1})), _uZstack.getStackSizes(), *(hh.get(Symmetry{1, 0, 0, 0})),
          ypart)
 {
   nvs[0] = ypart.nys();
diff --git a/dynare++/kord/korder_stoch.hh b/dynare++/kord/korder_stoch.hh
index 56ac74cea..65705afe3 100644
--- a/dynare++/kord/korder_stoch.hh
+++ b/dynare++/kord/korder_stoch.hh
@@ -87,7 +87,7 @@ IntegDerivs<t>::IntegDerivs(int r, const IntSequence &nvs, const _Tgss &g, const
       for (int i = 0; i <= d; i++)
         {
           int p = d-i;
-          Symmetry sym(i, 0, 0, p);
+          Symmetry sym{i, 0, 0, p};
           _Ttensor *ten = new _Ttensor(r, TensorDimens(sym, nvs));
 
           // calculate derivative $h_{y^i\sigma^p}$
@@ -105,14 +105,14 @@ IntegDerivs<t>::IntegDerivs(int r, const IntSequence &nvs, const _Tgss &g, const
               for (int m = 0; i+m+n+k <= maxd; m++, mfac *= m)
                 {
                   double mult = (pow(at_sigma, m)*povern)/mfac;
-                  Symmetry sym_mn(i, m+n, 0, k);
+                  Symmetry sym_mn{i, m+n, 0, k};
                   if (m+n == 0 && g.check(sym_mn))
                     ten->add(mult, *(g.get(sym_mn)));
                   if (m+n > 0 && KOrder::is_even(m+n) && g.check(sym_mn))
                     {
                       _Ttensor gtmp(*(g.get(sym_mn)));
                       gtmp.mult(mult);
-                      gtmp.contractAndAdd(1, *ten, *(mom.get(Symmetry(m+n))));
+                      gtmp.contractAndAdd(1, *ten, *(mom.get(Symmetry{m+n})));
                     }
                 }
             }
@@ -187,8 +187,8 @@ StochForwardDerivs<t>::StochForwardDerivs(const PartitionY &ypart, int nu,
       for (int i = 0; i <= d; i++)
         {
           int k = d-i;
-          if (g_int.check(Symmetry(i, 0, 0, k)))
-            ten->addSubTensor(*(g_int.get(Symmetry(i, 0, 0, k))));
+          if (g_int.check(Symmetry{i, 0, 0, k}))
+            ten->addSubTensor(*(g_int.get(Symmetry{i, 0, 0, k})));
         }
       g_int_sym.insert(ten);
     }
@@ -224,13 +224,13 @@ StochForwardDerivs<t>::StochForwardDerivs(const PartitionY &ypart, int nu,
   true_nvs[1] = nu; true_nvs[2] = nu;
   for (int d = 1; d <= maxd; d++)
     {
-      if (g_int_cent.check(Symmetry(d)))
+      if (g_int_cent.check(Symmetry{d}))
         {
           for (int i = 0; i <= d; i++)
             {
-              Symmetry sym(i, 0, 0, d-i);
+              Symmetry sym{i, 0, 0, d-i};
               IntSequence coor(sym, pp);
-              _Ttensor *ten = new _Ttensor(*(g_int_cent.get(Symmetry(d))), ss, coor,
+              _Ttensor *ten = new _Ttensor(*(g_int_cent.get(Symmetry{d})), ss, coor,
                                            TensorDimens(sym, true_nvs));
               this->insert(ten);
             }
@@ -274,7 +274,7 @@ GXContainer<_Ttype>::getType(int i, const Symmetry &s) const
   if (i == 2)
     return _Stype::zero;
   if (i == 3)
-    if (s == Symmetry(0, 0, 0, 1))
+    if (s == Symmetry{0, 0, 0, 1})
       return _Stype::unit;
     else
       return _Stype::zero;
@@ -319,12 +319,12 @@ ZXContainer<_Ttype>::getType(int i, const Symmetry &s) const
     else
       return _Stype::matrix;
   if (i == 2)
-    if (s == Symmetry(1, 0, 0, 0))
+    if (s == Symmetry{1, 0, 0, 0})
       return _Stype::unit;
     else
       return _Stype::zero;
   if (i == 3)
-    if (s == Symmetry(0, 1, 0, 0))
+    if (s == Symmetry{0, 1, 0, 0})
       return _Stype::unit;
     else
       return _Stype::zero;
diff --git a/dynare++/kord/tests.cc b/dynare++/kord/tests.cc
index f721417ec..6329754ad 100644
--- a/dynare++/kord/tests.cc
+++ b/dynare++/kord/tests.cc
@@ -287,7 +287,7 @@ TestRunnable::korder_unfold_fold(int maxdim, int unfold_dim,
   for (int d = 1; d <= maxdim; d++)
     {
       printf("\ttensor fill for dim=%d is:   %3.2f %%\n",
-             d, c.get(Symmetry(d))->getFillFactor()*100.0);
+             d, c.get(Symmetry{d})->getFillFactor()*100.0);
     }
   Journal jr("out.txt");
   KOrder kord(nstat, npred, nboth, nforw, c, gy, gu, v, jr);
diff --git a/dynare++/src/dynare3.cc b/dynare++/src/dynare3.cc
index 86d3ee12f..279c31003 100644
--- a/dynare++/src/dynare3.cc
+++ b/dynare++/src/dynare3.cc
@@ -353,7 +353,7 @@ DynareDerEvalLoader::DynareDerEvalLoader(const ogp::FineAtoms &a,
 void
 DynareDerEvalLoader::load(int i, int iord, const int *vars, double res)
 {
-  FSSparseTensor *t = md.get(Symmetry(iord));
+  FSSparseTensor *t = md.get(Symmetry{iord});
   IntSequence s(iord, 0);
   for (int j = 0; j < iord; j++)
     s[j] = atoms.get_pos_of_all(vars[j]);
diff --git a/dynare++/tl/cc/fs_tensor.hh b/dynare++/tl/cc/fs_tensor.hh
index 34e07e232..5417841c2 100644
--- a/dynare++/tl/cc/fs_tensor.hh
+++ b/dynare++/tl/cc/fs_tensor.hh
@@ -73,7 +73,7 @@ public:
   Symmetry
   getSym() const
   {
-    return Symmetry(dimen());
+    return Symmetry{dimen()};
   }
 
   int getOffset(const IntSequence &v) const override;
@@ -117,7 +117,7 @@ public:
   Symmetry
   getSym() const
   {
-    return Symmetry(dimen());
+    return Symmetry{dimen()};
   }
 
   int getOffset(const IntSequence &v) const override;
diff --git a/dynare++/tl/cc/gs_tensor.cc b/dynare++/tl/cc/gs_tensor.cc
index 5e53c93ce..cec65663c 100644
--- a/dynare++/tl/cc/gs_tensor.cc
+++ b/dynare++/tl/cc/gs_tensor.cc
@@ -9,7 +9,7 @@
    |@<|TensorDimens| class declaration@>| for details. */
 TensorDimens::TensorDimens(const IntSequence &ss, const IntSequence &coor)
   : nvs(ss),
-    sym(ss.size(), ""),
+    sym(ss.size()),
     nvmax(coor.size(), 0)
 {
   TL_RAISE_IF(!coor.isSorted(),
diff --git a/dynare++/tl/cc/gs_tensor.hh b/dynare++/tl/cc/gs_tensor.hh
index d266bf58b..89bc51b45 100644
--- a/dynare++/tl/cc/gs_tensor.hh
+++ b/dynare++/tl/cc/gs_tensor.hh
@@ -64,7 +64,7 @@ public:
   {
   }
   TensorDimens(int nvar, int dimen)
-    : nvs(1), sym(dimen), nvmax(dimen, nvar)
+    : nvs(1), sym{dimen}, nvmax(dimen, nvar)
   {
     nvs[0] = nvar;
   }
diff --git a/dynare++/tl/cc/int_sequence.hh b/dynare++/tl/cc/int_sequence.hh
index 11467f135..58dfd7e64 100644
--- a/dynare++/tl/cc/int_sequence.hh
+++ b/dynare++/tl/cc/int_sequence.hh
@@ -27,10 +27,16 @@
 #include <vector>
 #include <memory>
 #include <algorithm>
+#include <initializer_list>
 
 /* The implementation of |IntSequence| is straightforward. It has a pointer
    |data|, an |offset| integer indicating the beginning of the data relatively
-   to the pointer and a |length| of the sequence. */
+   to the pointer and a |length| of the sequence.
+
+   WARNING: IntSequence(n) and IntSequence{n} are not the same. The former
+   initializes a sequence of length n, while the latter constructs a sequence
+   of a single element equal to n. This is similar to the behaviour of
+   std::vector. */
 
 class Symmetry;
 class IntSequence
@@ -39,42 +45,52 @@ class IntSequence
   int length;
   int offset{0};
 public:
-  /* We have a constructor allocating a given length of data, constructor
-     allocating and then initializing all members to a given number, a copy
-     constructor, a conversion from |vector<int>|, a subsequence
-     constructor, a constructor used for calculating implied symmetry from
-     a more general symmetry and one equivalence class (see |Symmetry|
-     class). Finally we have a constructor which unfolds a sequence with
-     respect to a given symmetry and constructor which inserts a given
-     number to the ordered sequence or given number to a given position. */
-
+  // Constructor allocating a given length of (uninitialized) data
   explicit IntSequence(int l)
     : data{new int[l], [](int *arr) { delete[] arr; }}, length{l}
   {
   }
+  // Constructor allocating and then initializing all members to a given number
   IntSequence(int l, int n)
     : data{new int[l], [](int *arr) { delete[] arr; }}, length{l}
   {
     std::fill_n(data.get(), length, n);
   }
+  /* Constructor using an initializer list (gives the contents of the
+     IntSequence, similarly to std::vector) */
+  IntSequence(std::initializer_list<int> init)
+    : data{new int[init.size()], [](int *arr) { delete[] arr; }},
+      length{static_cast<int>(init.size())}
+  {
+    std::copy(init.begin(), init.end(), data.get());
+  }
+  // Copy constructor
   IntSequence(const IntSequence &s)
     : data{new int[s.length], [](int *arr) { delete[] arr; }}, length{s.length}
   {
     std::copy_n(s.data.get()+s.offset, length, data.get());
   }
+  // Move constructor
   IntSequence(IntSequence &&s) = default;
+  // Subsequence constructor (which shares the data pointer)
   IntSequence(IntSequence &s, int i1, int i2)
     : data{s.data}, length{i2-i1}, offset{s.offset+i1}
   {
   }
+  // Subsequence constructor (without pointer sharing)
   IntSequence(const IntSequence &s, int i1, int i2)
     : data{new int[i2-i1], [](int *arr) { delete[] arr; }}, length{i2-i1}
   {
     std::copy_n(s.data.get()+s.offset+i1, length, data.get());
   }
+  /* Constructor used for calculating implied symmetry from a more general
+     symmetry and one equivalence class */
   IntSequence(const Symmetry &sy, const std::vector<int> &se);
+  // Unfolds a given integer sequence with respect to a given symmetry
   IntSequence(const Symmetry &sy, const IntSequence &se);
+  // Inserts an element in an ordered sequence
   IntSequence(int i, const IntSequence &s);
+  // Inserts an element at a given position
   IntSequence(int i, const IntSequence &s, int pos);
 
   const IntSequence &operator=(const IntSequence &s);
diff --git a/dynare++/tl/cc/rfs_tensor.hh b/dynare++/tl/cc/rfs_tensor.hh
index 26717a579..adfd013c0 100644
--- a/dynare++/tl/cc/rfs_tensor.hh
+++ b/dynare++/tl/cc/rfs_tensor.hh
@@ -73,7 +73,7 @@ public:
   Symmetry
   getSym() const
   {
-    return Symmetry(dimen());
+    return Symmetry{dimen()};
   }
 };
 
@@ -113,7 +113,7 @@ public:
   Symmetry
   getSym() const
   {
-    return Symmetry(dimen());
+    return Symmetry{dimen()};
   }
 };
 
diff --git a/dynare++/tl/cc/sparse_tensor.cc b/dynare++/tl/cc/sparse_tensor.cc
index 6e3a1d129..5c3d59b91 100644
--- a/dynare++/tl/cc/sparse_tensor.cc
+++ b/dynare++/tl/cc/sparse_tensor.cc
@@ -116,7 +116,7 @@ SparseTensor::print() const
 
 FSSparseTensor::FSSparseTensor(int d, int nvar, int r)
   : SparseTensor(d, r, FFSTensor::calcMaxOffset(nvar, d)),
-    nv(nvar), sym(d)
+    nv(nvar), sym{d}
 {
 }
 
diff --git a/dynare++/tl/cc/stack_container.hh b/dynare++/tl/cc/stack_container.hh
index 37606a082..b40b6774d 100644
--- a/dynare++/tl/cc/stack_container.hh
+++ b/dynare++/tl/cc/stack_container.hh
@@ -284,8 +284,8 @@ public:
   multAndAdd(int dim, const TensorContainer<FSSparseTensor> &c,
              FGSTensor &out) const
   {
-    if (c.check(Symmetry(dim)))
-      multAndAdd(*(c.get(Symmetry(dim))), out);
+    if (c.check(Symmetry{dim}))
+      multAndAdd(*(c.get(Symmetry{dim})), out);
   }
   void multAndAdd(const FSSparseTensor &t, FGSTensor &out) const;
   void multAndAdd(int dim, const FGSContainer &c, FGSTensor &out) const;
@@ -314,8 +314,8 @@ public:
   multAndAdd(int dim, const TensorContainer<FSSparseTensor> &c,
              UGSTensor &out) const
   {
-    if (c.check(Symmetry(dim)))
-      multAndAdd(*(c.get(Symmetry(dim))), out);
+    if (c.check(Symmetry{dim}))
+      multAndAdd(*(c.get(Symmetry{dim})), out);
   }
   void multAndAdd(const FSSparseTensor &t, UGSTensor &out) const;
   void multAndAdd(int dim, const UGSContainer &c, UGSTensor &out) const;
@@ -370,12 +370,12 @@ public:
       else
         return _Stype::matrix;
     if (i == 2)
-      if (s == Symmetry(1, 0, 0, 0))
+      if (s == Symmetry{1, 0, 0, 0})
         return _Stype::unit;
       else
         return _Stype::zero;
     if (i == 3)
-      if (s == Symmetry(0, 1, 0, 0))
+      if (s == Symmetry{0, 1, 0, 0})
         return _Stype::unit;
       else
         return _Stype::zero;
@@ -446,19 +446,19 @@ public:
   getType(int i, const Symmetry &s) const override
   {
     if (i == 0)
-      if (s[2] > 0 || s == Symmetry(0, 0, 0, 1))
+      if (s[2] > 0 || s == Symmetry{0, 0, 0, 1})
         return _Stype::zero;
       else
         return _Stype::matrix;
     if (i == 1)
-      if (s == Symmetry(0, 0, 1, 0))
+      if (s == Symmetry{0, 0, 1, 0})
         return _Stype::unit;
       else
         return _Stype::zero;
     if (i == 2)
       return _Stype::zero;
     if (i == 3)
-      if (s == Symmetry(0, 0, 0, 1))
+      if (s == Symmetry{0, 0, 0, 1})
         return _Stype::unit;
       else
         return _Stype::zero;
diff --git a/dynare++/tl/cc/symmetry.cc b/dynare++/tl/cc/symmetry.cc
index ac54b841b..9b94df492 100644
--- a/dynare++/tl/cc/symmetry.cc
+++ b/dynare++/tl/cc/symmetry.cc
@@ -3,7 +3,7 @@
 #include "symmetry.hh"
 #include "permutation.hh"
 
-#include <cstdio>
+#include <iostream>
 
 /* Construct symmetry as numbers of successively equal items in the sequence. */
 
@@ -57,28 +57,18 @@ Symmetry::isFull() const
    beginning as subordinal |symiterator|. */
 
 symiterator::symiterator(SymmetrySet &ss)
-  : s(ss), subit(nullptr), subs(nullptr), end_flag(false)
+  : s(ss), end_flag(false)
 {
   s.sym()[0] = 0;
   if (s.size() == 2)
-    {
-      s.sym()[1] = s.dimen();
-    }
+    s.sym()[1] = s.dimen();
   else
     {
-      subs = new SymmetrySet(s, s.dimen());
-      subit = new symiterator(*subs);
+      subs = std::make_unique<SymmetrySet>(s, s.dimen());
+      subit = std::make_unique<symiterator>(*subs);
     }
 }
 
-symiterator::~symiterator()
-{
-  if (subit)
-    delete subit;
-  if (subs)
-    delete subs;
-}
-
 /* Here we move to the next symmetry. We do so only, if we are not at
    the end. If length is 2, we increase lower index and decrease upper
    index, otherwise we increase the subordinal symmetry. If we got to the
@@ -101,11 +91,9 @@ symiterator::operator++()
           ++(*subit);
           if (subit->isEnd())
             {
-              delete subit;
-              delete subs;
               s.sym()[0]++;
-              subs = new SymmetrySet(s, s.dimen()-s.sym()[0]);
-              subit = new symiterator(*subs);
+              subs = std::make_unique<SymmetrySet>(s, s.dimen()-s.sym()[0]);
+              subit = std::make_unique<symiterator>(*subs);
             }
         }
       if (s.sym()[0] == s.dimen()+1)
@@ -117,9 +105,7 @@ symiterator::operator++()
 InducedSymmetries::InducedSymmetries(const Equivalence &e, const Symmetry &s)
 {
   for (const auto & i : e)
-    {
-      push_back(Symmetry(s, i));
-    }
+    emplace_back(s, i);
 }
 
 // |InducedSymmetries| permuted constructor code
@@ -129,7 +115,7 @@ InducedSymmetries::InducedSymmetries(const Equivalence &e, const Permutation &p,
   for (int i = 0; i < e.numClasses(); i++)
     {
       auto it = e.find(p.getMap()[i]);
-      push_back(Symmetry(s, *it));
+      emplace_back(s, *it);
     }
 }
 
@@ -138,7 +124,7 @@ InducedSymmetries::InducedSymmetries(const Equivalence &e, const Permutation &p,
 void
 InducedSymmetries::print() const
 {
-  printf("Induced symmetries: %lu\n", (unsigned long) size());
+  std::cout << "Induced symmetries: " << size() << std::endl;
   for (unsigned int i = 0; i < size(); i++)
     operator[](i).print();
 }
diff --git a/dynare++/tl/cc/symmetry.hh b/dynare++/tl/cc/symmetry.hh
index eff2d84d0..2a4976ce7 100644
--- a/dynare++/tl/cc/symmetry.hh
+++ b/dynare++/tl/cc/symmetry.hh
@@ -46,63 +46,49 @@
 
 #include <list>
 #include <vector>
+#include <initializer_list>
+#include <utility>
+#include <memory>
 
 /* Clear. The method |isFull| returns true if and only if the symmetry
-   allows for any permutation of indices. */
+   allows for any permutation of indices.
+
+   WARNING: Symmetry(n) and Symmetry{n} are not the same. The former
+   initializes a symmetry of n elements, while the latter is a full symmetry of
+   order n. This is similar to the behaviour of std::vector. */
 
 class Symmetry : public IntSequence
 {
 public:
-  /* We provide three constructors for symmetries of the form $y^n$,
-     $y^nu^m$, $y^nu^m\sigma^k$. Also a copy constructor, and finally a
-     constructor of implied symmetry for a symmetry and an equivalence
-     class. It is already implemented in |IntSequence| so we only call
-     appropriate constructor of |IntSequence|. We also provide the
-     subsymmetry, which takes the given length of symmetry from the end.
-
-     The last constructor constructs a symmetry from an integer sequence
-     (supposed to be ordered) as a symmetry counting successively equal
-     items. For instance the sequence $(a,a,a,b,c,c,d,d,d,d)$ produces
-     symmetry $(3,1,2,4)$. */
-  Symmetry(int len, const char *dummy)
+  // Constructor allocating a given length of (zero-initialized) data
+  Symmetry(int len)
     : IntSequence(len, 0)
   {
   }
-  Symmetry(int i1)
-    : IntSequence(1, i1)
+  /* Constructor using an initializer list, that gives the contents of the
+     Symmetry. Used for symmetries of the form $y^n$, $y^n u^m$, $y^nu^m\sigma^k$ */
+  Symmetry(std::initializer_list<int> init)
+    : IntSequence(std::move(init))
   {
   }
-  Symmetry(int i1, int i2)
-    : IntSequence(2)
-  {
-    operator[](0) = i1;
-    operator[](1) = i2;
-  }
-  Symmetry(int i1, int i2, int i3)
-    : IntSequence(3)
-  {
-    operator[](0) = i1;
-    operator[](1) = i2;
-    operator[](2) = i3;
-  }
-  Symmetry(int i1, int i2, int i3, int i4)
-    : IntSequence(4)
-  {
-    operator[](0) = i1;
-    operator[](1) = i2;
-    operator[](2) = i3;
-    operator[](3) = i4;
-  }
+  // Copy constructor
   Symmetry(const Symmetry &s) = default;
+  // Move constructor
   Symmetry(Symmetry &&s) = default;
+  // Constructor of implied symmetry for a symmetry and an equivalence class
   Symmetry(const Symmetry &s, const OrdSequence &cl)
     : IntSequence(s, cl.getData())
   {
   }
+  /* Subsymmetry, which takes the given length of symmetry from the end (shares
+     data pointer) */
   Symmetry(Symmetry &s, int len)
     : IntSequence(s, s.size()-len, s.size())
   {
   }
+  /* Constructs a symmetry from an integer sequence (supposed to be ordered) as
+     a symmetry counting successively equal items. For instance the sequence
+     $(a,a,a,b,c,c,d,d,d,d)$ produces symmetry $(3,1,2,4)$. */
   Symmetry(const IntSequence &s);
 
   int
@@ -150,7 +136,7 @@ class SymmetrySet
   int dim;
 public:
   SymmetrySet(int d, int length)
-    : run(length, ""), dim(d)
+    : run(length), dim(d)
   {
   }
   SymmetrySet(SymmetrySet &s, int d)
@@ -184,7 +170,7 @@ public:
    to the |SymmetrySet| only to know dimension and for access of its
    symmetry storage. Further we have pointers to subordinal |symiterator|
    and its |SymmetrySet|. These are pointers, since the recursion ends at
-   length equal to 2, in which case these pointers are |NULL|.
+   length equal to 2, in which case these pointers are uninitialized.
 
    The constructor creates the iterator which initializes to the first
    symmetry (beginning). */
@@ -192,12 +178,12 @@ public:
 class symiterator
 {
   SymmetrySet &s;
-  symiterator *subit;
-  SymmetrySet *subs;
+  std::unique_ptr<symiterator> subit;
+  std::unique_ptr<SymmetrySet> subs;
   bool end_flag;
 public:
   symiterator(SymmetrySet &ss);
-  ~symiterator();
+  ~symiterator() = default;
   symiterator &operator++();
   bool
   isEnd() const
diff --git a/dynare++/tl/cc/t_polynomial.hh b/dynare++/tl/cc/t_polynomial.hh
index f1d258ee2..4130459be 100644
--- a/dynare++/tl/cc/t_polynomial.hh
+++ b/dynare++/tl/cc/t_polynomial.hh
@@ -168,14 +168,14 @@ public:
         const _Stype &xpow = pwp.getNext((const _Stype *) nullptr);
         for (int j = 0; j <= tp.maxdim-i; j++)
           {
-            if (tp.check(Symmetry(i+j)))
+            if (tp.check(Symmetry{i+j}))
               {
                 // initialize |ten| of dimension |j|
                 /* The pointer |ten| is either a new tensor or got from |this| container. */
                 _Ttype *ten;
-                if (_Tparent::check(Symmetry(j)))
+                if (_Tparent::check(Symmetry{j}))
                   {
-                    ten = _Tparent::get(Symmetry(j));
+                    ten = _Tparent::get(Symmetry{j});
                   }
                 else
                   {
@@ -184,9 +184,9 @@ public:
                     insert(ten);
                   }
 
-                Symmetry sym(i, j);
+                Symmetry sym{i, j};
                 IntSequence coor(sym, pp);
-                _TGStype slice(*(tp.get(Symmetry(i+j))), ss, coor, TensorDimens(sym, ss));
+                _TGStype slice(*(tp.get(Symmetry{i+j})), ss, coor, TensorDimens(sym, ss));
                 slice.mult(Tensor::noverk(i+j, j));
                 _TGStype tmp(*ten);
                 slice.contractAndAdd(0, tmp, xpow);
@@ -200,15 +200,15 @@ public:
        simple addition. */
     for (int j = 0; j <= tp.maxdim; j++)
       {
-        if (tp.check(Symmetry(j)))
+        if (tp.check(Symmetry{j}))
           {
 
             // initialize |ten| of dimension |j|
             /* Same code as above */
             _Ttype *ten;
-            if (_Tparent::check(Symmetry(j)))
+            if (_Tparent::check(Symmetry{j}))
               {
-                ten = _Tparent::get(Symmetry(j));
+                ten = _Tparent::get(Symmetry{j});
               }
             else
               {
@@ -217,9 +217,9 @@ public:
                 insert(ten);
               }
 
-            Symmetry sym(0, j);
+            Symmetry sym{0, j};
             IntSequence coor(sym, pp);
-            _TGStype slice(*(tp.get(Symmetry(j))), ss, coor, TensorDimens(sym, ss));
+            _TGStype slice(*(tp.get(Symmetry{j})), ss, coor, TensorDimens(sym, ss));
             ten->add(1.0, slice);
           }
       }
@@ -247,8 +247,8 @@ public:
   void
   evalTrad(Vector &out, const ConstVector &v) const
   {
-    if (_Tparent::check(Symmetry(0)))
-      out = _Tparent::get(Symmetry(0))->getData();
+    if (_Tparent::check(Symmetry{0}))
+      out = _Tparent::get(Symmetry{0})->getData();
     else
       out.zeros();
 
@@ -256,7 +256,7 @@ public:
     for (int d = 1; d <= maxdim; d++)
       {
         const _Stype &p = pp.getNext((const _Stype *) nullptr);
-        Symmetry cs(d);
+        Symmetry cs{d};
         if (_Tparent::check(cs))
           {
             const _Ttype *t = _Tparent::get(cs);
@@ -271,8 +271,8 @@ public:
   void
   evalHorner(Vector &out, const ConstVector &v) const
   {
-    if (_Tparent::check(Symmetry(0)))
-      out = _Tparent::get(Symmetry(0))->getData();
+    if (_Tparent::check(Symmetry{0}))
+      out = _Tparent::get(Symmetry{0})->getData();
     else
       out.zeros();
 
@@ -281,12 +281,12 @@ public:
 
     _Ttype *last;
     if (maxdim == 1)
-      last = new _Ttype(*(_Tparent::get(Symmetry(1))));
+      last = new _Ttype(*(_Tparent::get(Symmetry{1})));
     else
-      last = new _Ttype(*(_Tparent::get(Symmetry(maxdim))), v);
+      last = new _Ttype(*(_Tparent::get(Symmetry{maxdim})), v);
     for (int d = maxdim-1; d >= 1; d--)
       {
-        Symmetry cs(d);
+        Symmetry cs{d};
         if (_Tparent::check(cs))
           {
             const _Ttype *nt = _Tparent::get(cs);
@@ -337,9 +337,9 @@ public:
   {
     for (int d = 1; d <= maxdim; d++)
       {
-        if (_Tparent::check(Symmetry(d)))
+        if (_Tparent::check(Symmetry{d}))
           {
-            _Ttype *ten = _Tparent::get(Symmetry(d));
+            _Ttype *ten = _Tparent::get(Symmetry{d});
             ten->mult((double) std::max((d-k), 0));
           }
       }
@@ -367,14 +367,14 @@ public:
     auto *res = new _Ttype(nrows(), nvars(), s);
     res->zeros();
 
-    if (_Tparent::check(Symmetry(s)))
-      res->add(1.0, *(_Tparent::get(Symmetry(s))));
+    if (_Tparent::check(Symmetry{s}))
+      res->add(1.0, *(_Tparent::get(Symmetry{s})));
 
     for (int d = s+1; d <= maxdim; d++)
       {
-        if (_Tparent::check(Symmetry(d)))
+        if (_Tparent::check(Symmetry{d}))
           {
-            const _Ttype &ltmp = *(_Tparent::get(Symmetry(d)));
+            const _Ttype &ltmp = *(_Tparent::get(Symmetry{d}));
             auto *last = new _Ttype(ltmp);
             for (int j = 0; j < d - s; j++)
               {
@@ -474,12 +474,12 @@ public:
     for (Tensor::index i = dum.begin(); i != dum.end(); ++i)
       {
         int d = i.getCoor().sum();
-        Symmetry symrun(_Ttype::dimen()-d, d);
+        Symmetry symrun{_Ttype::dimen()-d, d};
         _TGStype dumgs(0, TensorDimens(symrun, dumnvs));
-        if (pol.check(Symmetry(d)))
+        if (pol.check(Symmetry{d}))
           {
             TwoDMatrix subt(*this, offset, dumgs.ncols());
-            subt.add(1.0, *(pol.get(Symmetry(d))));
+            subt.add(1.0, *(pol.get(Symmetry{d})));
           }
         offset += dumgs.ncols();
       }
diff --git a/dynare++/tl/testing/factory.hh b/dynare++/tl/testing/factory.hh
index 97281eb2a..f7b78c21a 100644
--- a/dynare++/tl/testing/factory.hh
+++ b/dynare++/tl/testing/factory.hh
@@ -51,7 +51,7 @@ public:
         if (symnum == 1)
           {
             // full symmetry
-            Symmetry sym(dim);
+            Symmetry sym{dim};
             auto *t = make<_Ttype>(r, sym, nvs);
             res->insert(t);
           }
@@ -60,7 +60,7 @@ public:
             // general symmetry
             for (int i = 0; i <= dim; i++)
               {
-                Symmetry sym(i, dim-i);
+                Symmetry sym{i, dim-i};
                 auto *t = make<_Ttype>(r, sym, nvs);
                 res->insert(t);
               }
diff --git a/dynare++/tl/testing/monoms.cc b/dynare++/tl/testing/monoms.cc
index 3fbe378c0..391f146d0 100644
--- a/dynare++/tl/testing/monoms.cc
+++ b/dynare++/tl/testing/monoms.cc
@@ -121,7 +121,7 @@ FGSTensor *
 Monom1Vector::deriv(int dim) const
 {
   FGSTensor *res
-    = new FGSTensor(len, TensorDimens(Symmetry(dim), IntSequence(1, nx)));
+    = new FGSTensor(len, TensorDimens(Symmetry{dim}, IntSequence(1, nx)));
   for (Tensor::index it = res->begin(); it != res->end(); ++it)
     {
       Vector outcol{res->getCol(*it)};
@@ -229,7 +229,7 @@ Monom2Vector::deriv(int maxdim) const
       for (int ydim = 0; ydim <= dim; ydim++)
         {
           int udim = dim - ydim;
-          Symmetry s(ydim, udim);
+          Symmetry s{ydim, udim};
           res->insert(deriv(s));
         }
     }
@@ -409,7 +409,7 @@ Monom4Vector::deriv(int dim) const
   FFSTensor dummy(0, nx1+nx2+nx3+nx4, dim);
   for (Tensor::index run = dummy.begin(); run != dummy.end(); ++run)
     {
-      Symmetry ind_sym(0, 0, 0, 0);
+      Symmetry ind_sym{0, 0, 0, 0};
       IntSequence ind(run.getCoor());
       for (int i = 0; i < ind.size(); i++)
         {
diff --git a/dynare++/tl/testing/tests.cc b/dynare++/tl/testing/tests.cc
index 558ca4090..d574fcb01 100644
--- a/dynare++/tl/testing/tests.cc
+++ b/dynare++/tl/testing/tests.cc
@@ -226,7 +226,7 @@ TestRunnable::dense_prod(const Symmetry &bsym, const IntSequence &bnvs,
   FGSContainer *cont
     = f.makeCont<FGSTensor, FGSContainer>(hnv, bnvs, bsym.dimen()-hdim+1);
   auto *fh
-    = f.make<FGSTensor>(rows, Symmetry(hdim), IntSequence(1, hnv));
+    = f.make<FGSTensor>(rows, Symmetry{hdim}, IntSequence(1, hnv));
   UGSTensor uh(*fh);
   FGSTensor fb(rows, TensorDimens(bsym, bnvs));
   fb.getData().zeros();
@@ -275,7 +275,7 @@ TestRunnable::folded_monomial(int ng, int nx, int ny, int nu, int dim)
   double maxnorm = 0;
   for (int ydim = 0; ydim <= dim; ydim++)
     {
-      Symmetry s(ydim, dim-ydim);
+      Symmetry s{ydim, dim-ydim};
       printf("\tSymmetry: "); s.print();
       FGSTensor res(ng, TensorDimens(s, nvs));
       res.getData().zeros();
@@ -314,7 +314,7 @@ TestRunnable::unfolded_monomial(int ng, int nx, int ny, int nu, int dim)
   double maxnorm = 0;
   for (int ydim = 0; ydim <= dim; ydim++)
     {
-      Symmetry s(ydim, dim-ydim);
+      Symmetry s{ydim, dim-ydim};
       printf("\tSymmetry: "); s.print();
       UGSTensor res(ng, TensorDimens(s, nvs));
       res.getData().zeros();
@@ -602,7 +602,7 @@ public:
   bool
   run() const override
   {
-    Symmetry s(2, 3);
+    Symmetry s{2, 3};
     IntSequence nvs(2); nvs[0] = 4; nvs[1] = 2;
     return index_forward<FGSTensor>(s, nvs);
   }
@@ -618,7 +618,7 @@ public:
   bool
   run() const override
   {
-    Symmetry s(2, 3);
+    Symmetry s{2, 3};
     IntSequence nvs(2); nvs[0] = 4; nvs[1] = 2;
     return index_forward<UGSTensor>(s, nvs);
   }
@@ -634,7 +634,7 @@ public:
   bool
   run() const override
   {
-    Symmetry s(2, 3, 2);
+    Symmetry s{2, 3, 2};
     IntSequence nvs(3); nvs[0] = 5; nvs[1] = 2; nvs[2] = 2;
     return index_forward<FGSTensor>(s, nvs);
   }
@@ -650,7 +650,7 @@ public:
   bool
   run() const override
   {
-    Symmetry s(2, 3, 2);
+    Symmetry s{2, 3, 2};
     IntSequence nvs(3); nvs[0] = 5; nvs[1] = 2; nvs[2] = 2;
     return index_forward<UGSTensor>(s, nvs);
   }
@@ -666,7 +666,7 @@ public:
   bool
   run() const override
   {
-    Symmetry s(1, 1, 3);
+    Symmetry s{1, 1, 3};
     IntSequence nvs(3); nvs[0] = 3; nvs[1] = 3; nvs[2] = 2;
     return index_backward<FGSTensor>(s, nvs);
   }
@@ -682,7 +682,7 @@ public:
   bool
   run() const override
   {
-    Symmetry s(2, 3, 2);
+    Symmetry s{2, 3, 2};
     IntSequence nvs(3); nvs[0] = 4; nvs[1] = 2; nvs[2] = 4;
     return index_backward<FGSTensor>(s, nvs);
   }
@@ -698,7 +698,7 @@ public:
   bool
   run() const override
   {
-    Symmetry s(1, 1, 3);
+    Symmetry s{1, 1, 3};
     IntSequence nvs(3); nvs[0] = 3; nvs[1] = 3; nvs[2] = 2;
     return index_backward<UGSTensor>(s, nvs);
   }
@@ -714,7 +714,7 @@ public:
   bool
   run() const override
   {
-    Symmetry s(2, 3, 2);
+    Symmetry s{2, 3, 2};
     IntSequence nvs(3); nvs[0] = 4; nvs[1] = 2; nvs[2] = 4;
     return index_backward<UGSTensor>(s, nvs);
   }
@@ -730,7 +730,7 @@ public:
   bool
   run() const override
   {
-    Symmetry s(2, 3);
+    Symmetry s{2, 3};
     IntSequence nvs(2); nvs[0] = 4; nvs[1] = 2;
     return index_offset<FGSTensor>(s, nvs);
   }
@@ -746,7 +746,7 @@ public:
   bool
   run() const override
   {
-    Symmetry s(2, 3);
+    Symmetry s{2, 3};
     IntSequence nvs(2); nvs[0] = 4; nvs[1] = 2;
     return index_offset<UGSTensor>(s, nvs);
   }
@@ -762,7 +762,7 @@ public:
   bool
   run() const override
   {
-    Symmetry s(2, 3, 2);
+    Symmetry s{2, 3, 2};
     IntSequence nvs(3); nvs[0] = 5; nvs[1] = 2; nvs[2] = 2;
     return index_offset<FGSTensor>(s, nvs);
   }
@@ -778,7 +778,7 @@ public:
   bool
   run() const override
   {
-    Symmetry s(2, 3, 2);
+    Symmetry s{2, 3, 2};
     IntSequence nvs(3); nvs[0] = 5; nvs[1] = 2; nvs[2] = 2;
     return index_offset<UGSTensor>(s, nvs);
   }
@@ -808,7 +808,7 @@ public:
   bool
   run() const override
   {
-    Symmetry s(1, 2, 2);
+    Symmetry s{1, 2, 2};
     IntSequence nvs(3); nvs[0] = 3; nvs[1] = 3; nvs[2] = 2;
     return gs_fold_unfold(5, s, nvs);
   }
@@ -838,7 +838,7 @@ public:
   bool
   run() const override
   {
-    Symmetry s(2, 1, 2);
+    Symmetry s{2, 1, 2};
     IntSequence nvs(3); nvs[0] = 6; nvs[1] = 2; nvs[2] = 6;
     return gs_fold_unfold(5, s, nvs);
   }
@@ -883,7 +883,7 @@ public:
   run() const override
   {
     IntSequence bnvs(2); bnvs[0] = 3; bnvs[1] = 2;
-    return dense_prod(Symmetry(1, 2), bnvs, 2, 3, 2);
+    return dense_prod(Symmetry{1, 2}, bnvs, 2, 3, 2);
   }
 };
 
@@ -898,7 +898,7 @@ public:
   run() const override
   {
     IntSequence bnvs(2); bnvs[0] = 10; bnvs[1] = 7;
-    return dense_prod(Symmetry(2, 3), bnvs, 3, 15, 10);
+    return dense_prod(Symmetry{2, 3}, bnvs, 3, 15, 10);
   }
 };
 
@@ -913,7 +913,7 @@ public:
   run() const override
   {
     IntSequence bnvs(2); bnvs[0] = 13; bnvs[1] = 11;
-    return dense_prod(Symmetry(3, 2), bnvs, 3, 20, 20);
+    return dense_prod(Symmetry{3, 2}, bnvs, 3, 20, 20);
   }
 };
 
diff --git a/mex/sources/k_order_perturbation/k_ord_dynare.cc b/mex/sources/k_order_perturbation/k_ord_dynare.cc
index 36d2b2634..de3389239 100644
--- a/mex/sources/k_order_perturbation/k_ord_dynare.cc
+++ b/mex/sources/k_order_perturbation/k_ord_dynare.cc
@@ -248,7 +248,7 @@ KordpDynare::populateDerivativesContainer(const TwoDMatrix &g, int ord, const st
     }
 
   // md container
-  md.remove(Symmetry(ord));
+  md.remove(Symmetry{ord});
   md.insert(mdTi);
   // No need to delete mdTi, it will be deleted by TensorContainer destructor
 }
diff --git a/mex/sources/k_order_perturbation/k_order_perturbation.cc b/mex/sources/k_order_perturbation/k_order_perturbation.cc
index f374e5641..f2581c79b 100644
--- a/mex/sources/k_order_perturbation/k_order_perturbation.cc
+++ b/mex/sources/k_order_perturbation/k_order_perturbation.cc
@@ -294,18 +294,18 @@ extern "C" {
                 for (int i = 0; i < 12; ++i)
                   c_fieldnames[i] = fieldnames[i].c_str();
                 plhs[ii] = mxCreateStructMatrix(1, 1, 12, c_fieldnames);
-                copy_derivatives(plhs[ii], Symmetry(1, 0, 0, 0), derivs, "gy");
-                copy_derivatives(plhs[ii], Symmetry(0, 1, 0, 0), derivs, "gu");
-                copy_derivatives(plhs[ii], Symmetry(2, 0, 0, 0), derivs, "gyy");
-                copy_derivatives(plhs[ii], Symmetry(0, 2, 0, 0), derivs, "guu");
-                copy_derivatives(plhs[ii], Symmetry(1, 1, 0, 0), derivs, "gyu");
-                copy_derivatives(plhs[ii], Symmetry(0, 0, 0, 2), derivs, "gss");
-                copy_derivatives(plhs[ii], Symmetry(3, 0, 0, 0), derivs, "gyyy");
-                copy_derivatives(plhs[ii], Symmetry(0, 3, 0, 0), derivs, "guuu");
-                copy_derivatives(plhs[ii], Symmetry(2, 1, 0, 0), derivs, "gyyu");
-                copy_derivatives(plhs[ii], Symmetry(1, 2, 0, 0), derivs, "gyuu");
-                copy_derivatives(plhs[ii], Symmetry(1, 0, 0, 2), derivs, "gyss");
-                copy_derivatives(plhs[ii], Symmetry(0, 1, 0, 2), derivs, "guss");
+                copy_derivatives(plhs[ii], Symmetry{1, 0, 0, 0}, derivs, "gy");
+                copy_derivatives(plhs[ii], Symmetry{0, 1, 0, 0}, derivs, "gu");
+                copy_derivatives(plhs[ii], Symmetry{2, 0, 0, 0}, derivs, "gyy");
+                copy_derivatives(plhs[ii], Symmetry{0, 2, 0, 0}, derivs, "guu");
+                copy_derivatives(plhs[ii], Symmetry{1, 1, 0, 0}, derivs, "gyu");
+                copy_derivatives(plhs[ii], Symmetry{0, 0, 0, 2}, derivs, "gss");
+                copy_derivatives(plhs[ii], Symmetry{3, 0, 0, 0}, derivs, "gyyy");
+                copy_derivatives(plhs[ii], Symmetry{0, 3, 0, 0}, derivs, "guuu");
+                copy_derivatives(plhs[ii], Symmetry{2, 1, 0, 0}, derivs, "gyyu");
+                copy_derivatives(plhs[ii], Symmetry{1, 2, 0, 0}, derivs, "gyuu");
+                copy_derivatives(plhs[ii], Symmetry{1, 0, 0, 2}, derivs, "gyss");
+                copy_derivatives(plhs[ii], Symmetry{0, 1, 0, 2}, derivs, "guss");
               }
           }
       }