diff --git a/dynare++/extern/matlab/dynare_simul.cc b/dynare++/extern/matlab/dynare_simul.cc
index 0bc2c70ff..f057ce738 100644
--- a/dynare++/extern/matlab/dynare_simul.cc
+++ b/dynare++/extern/matlab/dynare_simul.cc
@@ -102,8 +102,7 @@ extern "C" {
             ft.zeros();
             ConstTwoDMatrix gk_mat(ft.nrows(), ft.ncols(), ConstVector{gk});
             ft.add(1.0, gk_mat);
-            auto *ut = new UFSTensor(ft);
-            pol.insert(ut);
+            pol.insert(std::make_unique<UFSTensor>(ft));
           }
         // form the decision rule
         UnfoldDecisionRule
diff --git a/dynare++/kord/approximation.cc b/dynare++/kord/approximation.cc
index 3ed996d4b..fbda76025 100644
--- a/dynare++/kord/approximation.cc
+++ b/dynare++/kord/approximation.cc
@@ -1,5 +1,7 @@
 // Copyright 2005, Ondra Kamenik
 
+#include <utility>
+
 #include "kord_exception.hh"
 #include "approximation.hh"
 #include "first_order.hh"
@@ -261,8 +263,8 @@ Approximation::saveRuleDerivs(const FGSContainer &g)
   rule_ders_ss = new FGSContainer(4);
   for (auto & run : (*rule_ders))
     {
-      auto *ten = new FGSTensor(ypart.nstat+ypart.npred, ypart.nyss(), *(run.second));
-      rule_ders_ss->insert(ten);
+      auto ten = std::make_unique<FGSTensor>(ypart.nstat+ypart.npred, ypart.nyss(), *(run.second));
+      rule_ders_ss->insert(std::move(ten));
     }
 }
 
diff --git a/dynare++/kord/decision_rule.cc b/dynare++/kord/decision_rule.cc
index 930bf801c..8e9f0858b 100644
--- a/dynare++/kord/decision_rule.cc
+++ b/dynare++/kord/decision_rule.cc
@@ -35,9 +35,7 @@ FoldDecisionRule::FoldDecisionRule(const UnfoldDecisionRule &udr)
                                    udr.ypart, udr.nu, udr.ysteady)
 {
   for (const auto & it : udr)
-    {
-      insert(new ctraits<KOrder::fold>::Ttensym(*(it.second)));
-    }
+    insert(std::make_unique<ctraits<KOrder::fold>::Ttensym>(*(it.second)));
 }
 
 // |UnfoldDecisionRule| conversion from |FoldDecisionRule|
@@ -46,9 +44,7 @@ UnfoldDecisionRule::UnfoldDecisionRule(const FoldDecisionRule &fdr)
                                      fdr.ypart, fdr.nu, fdr.ysteady)
 {
   for (const auto & it : fdr)
-    {
-      insert(new ctraits<KOrder::unfold>::Ttensym(*(it.second)));
-    }
+    insert(std::make_unique<ctraits<KOrder::unfold>::Ttensym>(*(it.second)));
 }
 
 SimResults::~SimResults()
diff --git a/dynare++/kord/decision_rule.hh b/dynare++/kord/decision_rule.hh
index 749f7e175..dcc239c6a 100644
--- a/dynare++/kord/decision_rule.hh
+++ b/dynare++/kord/decision_rule.hh
@@ -198,7 +198,7 @@ DecisionRuleImpl<t>::fillTensors(const _Tg &g, double sigma)
   int dfact = 1;
   for (int d = 0; d <= g.getMaxDim(); d++, dfact *= d)
     {
-      auto *g_yud = new _Ttensym(ypart.ny(), ypart.nys()+nu, d);
+      auto g_yud = std::make_unique<_Ttensym>(ypart.ny(), ypart.nys()+nu, d);
       g_yud->zeros();
 
       // fill tensor of |g_yud| of dimension |d|
@@ -231,7 +231,7 @@ DecisionRuleImpl<t>::fillTensors(const _Tg &g, double sigma)
           g_yud->addSubTensor(tmp);
         }
 
-      this->insert(g_yud);
+      this->insert(std::move(g_yud));
     }
 }
 
@@ -261,9 +261,9 @@ DecisionRuleImpl<t>::centralize(const DecisionRuleImpl &dr)
   for (int d = 1; d <= dr.getMaxDim(); d++, dfac *= d)
     {
       pol.derivative(d-1);
-      _Ttensym *der = pol.evalPartially(d, dstate);
+      auto der = pol.evalPartially(d, dstate);
       der->mult(1.0/dfac);
-      this->insert(der);
+      this->insert(std::move(der));
     }
 }
 
@@ -567,7 +567,7 @@ DRFixPoint<t>::fillTensors(const _Tg &g, double sigma)
   int dfact = 1;
   for (int d = 0; d <= g.getMaxDim(); d++, dfact *= d)
     {
-      auto *g_yd = new _Ttensym(ypart.ny(), ypart.nys(), d);
+      auto g_yd = std::make_unique<_Ttensym>(ypart.ny(), ypart.nys(), d);
       g_yd->zeros();
       int kfact = 1;
       for (int k = 0; d+k <= g.getMaxDim(); k++, kfact *= k)
@@ -579,7 +579,7 @@ DRFixPoint<t>::fillTensors(const _Tg &g, double sigma)
               g_yd->add(mult, *ten);
             }
         }
-      this->insert(g_yd);
+      this->insert(std::move(g_yd));
     }
 }
 
@@ -613,7 +613,7 @@ DRFixPoint<t>::solveNewton(Vector &y)
 
   do
     {
-      _Ttensym *jacob = bigfder->evalPartially(1, sol);
+      auto jacob = bigfder->evalPartially(1, sol);
       bigf->evalHorner(delta, sol);
       if (newton_iter_last == 0)
         flastnorm = delta.getNorm();
@@ -646,7 +646,6 @@ DRFixPoint<t>::solveNewton(Vector &y)
           sol.add(-urelax, delta);
           delta_finite = delta.isFinite();
         }
-      delete jacob;
       newton_iter_last++;
       converged = delta_finite && fnorm < tol;
       flastnorm = fnorm;
diff --git a/dynare++/kord/first_order.hh b/dynare++/kord/first_order.hh
index c396570dc..6eee40ced 100644
--- a/dynare++/kord/first_order.hh
+++ b/dynare++/kord/first_order.hh
@@ -72,12 +72,14 @@ 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));
-    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->zeros(); ten->add(1.0, fo.gu);
-    this->insert(ten);
+    auto ten = std::make_unique<_Ttensor>(fo.ypart.ny(), TensorDimens(Symmetry{1, 0, 0, 0}, nvs));
+    ten->zeros();
+    ten->add(1.0, fo.gy);
+    this->insert(std::move(ten));
+    ten = std::make_unique<_Ttensor>(fo.ypart.ny(), TensorDimens(Symmetry{0, 1, 0, 0}, nvs));
+    ten->zeros();
+    ten->add(1.0, fo.gu);
+    this->insert(std::move(ten));
   }
 };
 
diff --git a/dynare++/kord/global_check.cc b/dynare++/kord/global_check.cc
index 08fdca3bc..1d3b01b7e 100644
--- a/dynare++/kord/global_check.cc
+++ b/dynare++/kord/global_check.cc
@@ -101,9 +101,9 @@ ResidFunction::setYU(const ConstVector &ys, const ConstVector &xx)
     }
   else
     {
-      auto *ten = new FFSTensor(hss->nrows(), hss->nvars(), 0);
+      auto ten = std::make_unique<FFSTensor>(hss->nrows(), hss->nvars(), 0);
       ten->getData() = ysteady_ss;
-      hss->insert(ten);
+      hss->insert(std::move(ten));
     }
 
 }
diff --git a/dynare++/kord/korder.cc b/dynare++/kord/korder.cc
index 553f1eb69..3d9022d23 100644
--- a/dynare++/kord/korder.cc
+++ b/dynare++/kord/korder.cc
@@ -265,21 +265,21 @@ 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));
+  auto tgy = std::make_unique<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));
+  insertDerivative<unfold>(std::move(tgy));
+  auto tgu = std::make_unique<UGSTensor>(ny, TensorDimens(Symmetry{0, 1, 0, 0}, nvs));
   tgu->getData() = gu.getData();
-  insertDerivative<unfold>(tgu);
+  insertDerivative<unfold>(std::move(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});
-  G<unfold>().insert(tGy);
-  UGSTensor *tGu = faaDiBrunoG<unfold>(Symmetry{0, 1, 0, 0});
-  G<unfold>().insert(tGu);
-  UGSTensor *tGup = faaDiBrunoG<unfold>(Symmetry{0, 0, 1, 0});
-  G<unfold>().insert(tGup);
+  auto tGy = faaDiBrunoG<unfold>(Symmetry{1, 0, 0, 0});
+  G<unfold>().insert(std::move(tGy));
+  auto tGu = faaDiBrunoG<unfold>(Symmetry{0, 1, 0, 0});
+  G<unfold>().insert(std::move(tGu));
+  auto tGup = faaDiBrunoG<unfold>(Symmetry{0, 0, 1, 0});
+  G<unfold>().insert(std::move(tGup));
 }
 
 // |KOrder::sylvesterSolve| unfolded specialization
@@ -346,8 +346,8 @@ KOrder::switchToFolded()
       {
         if (si[2] == 0 && g<unfold>().check(si))
           {
-            auto *ft = new FGSTensor(*(g<unfold>().get(si)));
-            insertDerivative<fold>(ft);
+            auto ft = std::make_unique<FGSTensor>(*(g<unfold>().get(si)));
+            insertDerivative<fold>(std::move(ft));
             if (dim > 1)
               {
                 gss<unfold>().remove(si);
@@ -357,8 +357,8 @@ KOrder::switchToFolded()
           }
         if (G<unfold>().check(si))
           {
-            auto *ft = new FGSTensor(*(G<unfold>().get(si)));
-            G<fold>().insert(ft);
+            auto ft = std::make_unique<FGSTensor>(*(G<unfold>().get(si)));
+            G<fold>().insert(std::move(ft));
             if (dim > 1)
               {
                 G<fold>().remove(si);
diff --git a/dynare++/kord/korder.hh b/dynare++/kord/korder.hh
index 5e3df787c..d4b51dcfa 100644
--- a/dynare++/kord/korder.hh
+++ b/dynare++/kord/korder.hh
@@ -361,14 +361,14 @@ public:
   }
 protected:
   template <int t>
-  void insertDerivative(_Ttensor *der);
+  void insertDerivative(std::unique_ptr<_Ttensor> der);
   template<int t>
   void sylvesterSolve(_Ttensor &der) const;
 
   template <int t>
-  _Ttensor *faaDiBrunoZ(const Symmetry &sym) const;
+  std::unique_ptr<_Ttensor> faaDiBrunoZ(const Symmetry &sym) const;
   template <int t>
-  _Ttensor *faaDiBrunoG(const Symmetry &sym) const;
+  std::unique_ptr<_Ttensor> faaDiBrunoG(const Symmetry &sym) const;
 
   template<int t>
   void recover_y(int i);
@@ -403,12 +403,13 @@ protected:
 
 template <int t>
 void
-KOrder::insertDerivative(_Ttensor *der)
+KOrder::insertDerivative(std::unique_ptr<_Ttensor> der)
 {
-  g<t>().insert(der);
-  gs<t>().insert(new _Ttensor(ypart.nstat, ypart.nys(), *der));
-  gss<t>().insert(new _Ttensor(ypart.nstat+ypart.npred,
-                               ypart.nyss(), *der));
+  auto der_ptr = der.get();
+  g<t>().insert(std::move(der));
+  gs<t>().insert(std::make_unique<_Ttensor>(ypart.nstat, ypart.nys(), *der_ptr));
+  gss<t>().insert(std::make_unique<_Ttensor>(ypart.nstat+ypart.npred,
+                                             ypart.nyss(), *der_ptr));
 }
 
 /* Here we implement Faa Di Bruno formula
@@ -419,12 +420,12 @@ KOrder::insertDerivative(_Ttensor *der)
    symmetry. */
 
 template <int t>
-_Ttensor *
+std::unique_ptr<_Ttensor>
 KOrder::faaDiBrunoZ(const Symmetry &sym) const
 {
   JournalRecordPair pa(journal);
   pa << "Faa Di Bruno Z container for " << sym << endrec;
-  _Ttensor *res = new _Ttensor(ny, TensorDimens(sym, nvs));
+  auto res = std::make_unique<_Ttensor>(ny, TensorDimens(sym, nvs));
   FaaDiBruno bruno(journal);
   bruno.calculate(Zstack<t>(), f, *res);
   return res;
@@ -434,13 +435,13 @@ KOrder::faaDiBrunoZ(const Symmetry &sym) const
    $g^{**}$ and $G$ stack. */
 
 template <int t>
-_Ttensor *
+std::unique_ptr<_Ttensor>
 KOrder::faaDiBrunoG(const Symmetry &sym) const
 {
   JournalRecordPair pa(journal);
   pa << "Faa Di Bruno G container for " << sym << endrec;
   TensorDimens tdims(sym, nvs);
-  auto *res = new _Ttensor(ypart.nyss(), tdims);
+  auto res = std::make_unique<_Ttensor>(ypart.nyss(), tdims);
   FaaDiBruno bruno(journal);
   bruno.calculate(Gstack<t>(), gss<t>(), *res);
   return res;
@@ -466,20 +467,21 @@ KOrder::recover_y(int i)
   JournalRecordPair pa(journal);
   pa << "Recovering symmetry " << sym << endrec;
 
-  _Ttensor *G_yi = faaDiBrunoG<t>(sym);
-  G<t>().insert(G_yi);
+  auto G_yi = faaDiBrunoG<t>(sym);
+  auto G_yi_ptr = G_yi.get();
+  G<t>().insert(std::move(G_yi));
 
-  _Ttensor *g_yi = faaDiBrunoZ<t>(sym);
+  auto g_yi = faaDiBrunoZ<t>(sym);
   g_yi->mult(-1.0);
 
   sylvesterSolve<t>(*g_yi);
 
-  insertDerivative<t>(g_yi);
+  insertDerivative<t>(std::move(g_yi));
 
   _Ttensor *gss_y = gss<t>().get(Symmetry{1, 0, 0, 0});
-  gs<t>().multAndAdd(*gss_y, *G_yi);
+  gs<t>().multAndAdd(*gss_y, *G_yi_ptr);
   _Ttensor *gss_yi = gss<t>().get(sym);
-  gs<t>().multAndAdd(*gss_yi, *G_yi);
+  gs<t>().multAndAdd(*gss_yi, *G_yi_ptr);
 }
 
 /* Here we solve $\left[F_{y^iu^j}\right]=0$ to obtain $g_{y^iu^j}$ for
@@ -501,15 +503,16 @@ KOrder::recover_yu(int i, int j)
   JournalRecordPair pa(journal);
   pa << "Recovering symmetry " << sym << endrec;
 
-  _Ttensor *G_yiuj = faaDiBrunoG<t>(sym);
-  G<t>().insert(G_yiuj);
+  auto G_yiuj = faaDiBrunoG<t>(sym);
+  auto G_yiuj_ptr = G_yiuj.get();
+  G<t>().insert(std::move(G_yiuj));
 
-  _Ttensor *g_yiuj = faaDiBrunoZ<t>(sym);
+  auto g_yiuj = faaDiBrunoZ<t>(sym);
   g_yiuj->mult(-1.0);
   matA.multInv(*g_yiuj);
-  insertDerivative<t>(g_yiuj);
+  insertDerivative<t>(std::move(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_ptr);
 }
 
 /* Here we solve
@@ -544,10 +547,11 @@ KOrder::recover_ys(int i, int j)
 
   if (is_even(j))
     {
-      _Ttensor *G_yisj = faaDiBrunoG<t>(sym);
-      G<t>().insert(G_yisj);
+      auto G_yisj = faaDiBrunoG<t>(sym);
+      auto G_yisj_ptr = G_yisj.get();
+      G<t>().insert(std::move(G_yisj));
 
-      _Ttensor *g_yisj = faaDiBrunoZ<t>(sym);
+      auto g_yisj = faaDiBrunoZ<t>(sym);
 
       {
         _Ttensor *D_ij = calcD_ik<t>(i, j);
@@ -566,10 +570,10 @@ KOrder::recover_ys(int i, int j)
 
       sylvesterSolve<t>(*g_yisj);
 
-      insertDerivative<t>(g_yisj);
+      insertDerivative<t>(std::move(g_yisj));
 
-      Gstack<t>().multAndAdd(1, gss<t>(), *G_yisj);
-      Gstack<t>().multAndAdd(i+j, gss<t>(), *G_yisj);
+      Gstack<t>().multAndAdd(1, gss<t>(), *G_yisj_ptr);
+      Gstack<t>().multAndAdd(i+j, gss<t>(), *G_yisj_ptr);
     }
 }
 
@@ -604,10 +608,11 @@ KOrder::recover_yus(int i, int j, int k)
 
   if (is_even(k))
     {
-      _Ttensor *G_yiujsk = faaDiBrunoG<t>(sym);
-      G<t>().insert(G_yiujsk);
+      auto G_yiujsk = faaDiBrunoG<t>(sym);
+      auto G_yiujsk_ptr = G_yiujsk.get();
+      G<t>().insert(std::move(G_yiujsk));
 
-      _Ttensor *g_yiujsk = faaDiBrunoZ<t>(sym);
+      auto g_yiujsk = faaDiBrunoZ<t>(sym);
 
       {
         _Ttensor *D_ijk = calcD_ijk<t>(i, j, k);
@@ -625,9 +630,9 @@ KOrder::recover_yus(int i, int j, int k)
       g_yiujsk->mult(-1.0);
 
       matA.multInv(*g_yiujsk);
-      insertDerivative<t>(g_yiujsk);
+      insertDerivative<t>(std::move(g_yiujsk));
 
-      Gstack<t>().multAndAdd(1, gss<t>(), *G_yiujsk);
+      Gstack<t>().multAndAdd(1, gss<t>(), *G_yiujsk_ptr);
     }
 }
 
@@ -671,10 +676,11 @@ KOrder::recover_s(int i)
 
   if (is_even(i))
     {
-      _Ttensor *G_si = faaDiBrunoG<t>(sym);
-      G<t>().insert(G_si);
+      auto G_si = faaDiBrunoG<t>(sym);
+      auto G_si_ptr = G_si.get();
+      G<t>().insert(std::move(G_si));
 
-      _Ttensor *g_si = faaDiBrunoZ<t>(sym);
+      auto g_si = faaDiBrunoZ<t>(sym);
 
       {
         _Ttensor *D_i = calcD_k<t>(i);
@@ -692,10 +698,10 @@ KOrder::recover_s(int i)
       g_si->mult(-1.0);
 
       matS.multInv(*g_si);
-      insertDerivative<t>(g_si);
+      insertDerivative<t>(std::move(g_si));
 
-      Gstack<t>().multAndAdd(1, gss<t>(), *G_si);
-      Gstack<t>().multAndAdd(i, gss<t>(), *G_si);
+      Gstack<t>().multAndAdd(1, gss<t>(), *G_si_ptr);
+      Gstack<t>().multAndAdd(i, gss<t>(), *G_si_ptr);
     }
 }
 
@@ -711,8 +717,8 @@ KOrder::fillG(int i, int j, int k)
     {
       if (is_even(k-m))
         {
-          _Ttensor *G_yiujupms = faaDiBrunoG<t>(Symmetry{i, j, m, k-m});
-          G<t>().insert(G_yiujupms);
+          auto G_yiujupms = faaDiBrunoG<t>(Symmetry{i, j, m, k-m});
+          G<t>().insert(std::move(G_yiujupms));
         }
     }
 }
@@ -732,9 +738,8 @@ KOrder::calcD_ijk(int i, int j, int k) const
   res->zeros();
   if (is_even(k))
     {
-      _Ttensor *tmp = faaDiBrunoZ<t>(Symmetry{i, j, k, 0});
+      auto tmp = faaDiBrunoZ<t>(Symmetry{i, j, k, 0});
       tmp->contractAndAdd(2, *res, *(m<t>().get(Symmetry{k})));
-      delete tmp;
     }
   return res;
 }
@@ -754,10 +759,9 @@ KOrder::calcE_ijk(int i, int j, int k) const
   res->zeros();
   for (int n = 2; n <= k-1; n += 2)
     {
-      _Ttensor *tmp = faaDiBrunoZ<t>(Symmetry{i, j, n, k-n});
+      auto tmp = faaDiBrunoZ<t>(Symmetry{i, j, n, k-n});
       tmp->mult((double) (PascalTriangle::noverk(k, n)));
       tmp->contractAndAdd(2, *res, *(m<t>().get(Symmetry{n})));
-      delete tmp;
     }
   return res;
 }
@@ -863,12 +867,11 @@ KOrder::check(int dim) const
   for (int i = 0; i <= dim; i++)
     {
       Symmetry sym{dim-i, i, 0,  0};
-      _Ttensor *r = faaDiBrunoZ<t>(sym);
+      auto r = faaDiBrunoZ<t>(sym);
       double err = r->getData().getMax();
       JournalRecord(journal) << "\terror for symmetry " << sym << "\tis " << err << endrec;
       if (err > maxerror)
         maxerror = err;
-      delete r;
     }
 
   // check for $F_{y^iu^ju'^k}+D_{ijk}+E_{ijk}=0$
@@ -880,7 +883,7 @@ KOrder::check(int dim) const
       if (i+j > 0 && k > 0)
         {
           Symmetry sym{i, j, 0, k};
-          _Ttensor *r = faaDiBrunoZ<t>(sym);
+          auto r = faaDiBrunoZ<t>(sym);
           _Ttensor *D_ijk = calcD_ijk<t>(i, j, k);
           r->add(1.0, *D_ijk);
           delete D_ijk;
@@ -889,12 +892,11 @@ KOrder::check(int dim) const
           delete E_ijk;
           double err = r->getData().getMax();
           JournalRecord(journal) << "\terror for symmetry " << sym << "\tis " << err << endrec;
-          delete r;
         }
     }
 
   // check for $F_{\sigma^i}+D_i+E_i=0
-  _Ttensor *r = faaDiBrunoZ<t>(Symmetry{0, 0, 0, dim});
+  auto r = faaDiBrunoZ<t>(Symmetry{0, 0, 0, dim});
   _Ttensor *D_k = calcD_k<t>(dim);
   r->add(1.0, *D_k);
   delete D_k;
@@ -906,7 +908,6 @@ KOrder::check(int dim) const
   JournalRecord(journal) << "\terror for symmetry " << sym << "\tis " << err << endrec;
   if (err > maxerror)
     maxerror = err;
-  delete r;
 
   return maxerror;
 }
diff --git a/dynare++/kord/korder_stoch.hh b/dynare++/kord/korder_stoch.hh
index 728a59b08..6f2bbfe56 100644
--- a/dynare++/kord/korder_stoch.hh
+++ b/dynare++/kord/korder_stoch.hh
@@ -22,6 +22,8 @@
    calculate $h$ as an extrapolation based on an approximation to $g$ at
    lower $\sigma$. */
 
+#include <memory>
+
 #include "korder.hh"
 #include "faa_di_bruno.hh"
 #include "journal.hh"
@@ -89,7 +91,7 @@ IntegDerivs<t>::IntegDerivs(int r, const IntSequence &nvs, const _Tgss &g, const
         {
           int p = d-i;
           Symmetry sym{i, 0, 0, p};
-          _Ttensor *ten = new _Ttensor(r, TensorDimens(sym, nvs));
+          auto ten = std::make_unique<_Ttensor>(r, TensorDimens(sym, nvs));
 
           // calculate derivative $h_{y^i\sigma^p}$
           /* This code calculates
@@ -118,7 +120,7 @@ IntegDerivs<t>::IntegDerivs(int r, const IntSequence &nvs, const _Tgss &g, const
                 }
             }
 
-          this->insert(ten);
+          this->insert(std::move(ten));
         }
     }
 }
@@ -183,7 +185,7 @@ StochForwardDerivs<t>::StochForwardDerivs(const PartitionY &ypart, int nu,
   _Tpol g_int_sym(r, ypart.nys()+1);
   for (int d = 1; d <= maxd; d++)
     {
-      auto *ten = new _Ttensym(r, ypart.nys()+1, d);
+      auto ten = std::make_unique<_Ttensym>(r, ypart.nys()+1, d);
       ten->zeros();
       for (int i = 0; i <= d; i++)
         {
@@ -191,7 +193,7 @@ StochForwardDerivs<t>::StochForwardDerivs(const PartitionY &ypart, int nu,
           if (g_int.check(Symmetry{i, 0, 0, k}))
             ten->addSubTensor(*(g_int.get(Symmetry{i, 0, 0, k})));
         }
-      g_int_sym.insert(ten);
+      g_int_sym.insert(std::move(ten));
     }
 
   // make |g_int_cent| the centralized polynomial about $(\bar y,\bar\sigma)$
@@ -210,8 +212,8 @@ StochForwardDerivs<t>::StochForwardDerivs(const PartitionY &ypart, int nu,
   for (int d = 1; d <= maxd; d++)
     {
       g_int_sym.derivative(d-1);
-      _Ttensym *der = g_int_sym.evalPartially(d, delta);
-      g_int_cent.insert(der);
+      auto der = g_int_sym.evalPartially(d, delta);
+      g_int_cent.insert(std::move(der));
     }
 
   // pull out general symmetry tensors from |g_int_cent|
@@ -231,9 +233,9 @@ StochForwardDerivs<t>::StochForwardDerivs(const PartitionY &ypart, int nu,
             {
               Symmetry sym{i, 0, 0, d-i};
               IntSequence coor(sym, pp);
-              _Ttensor *ten = new _Ttensor(*(g_int_cent.get(Symmetry{d})), ss, coor,
-                                           TensorDimens(sym, true_nvs));
-              this->insert(ten);
+              auto ten = std::make_unique<_Ttensor>(*(g_int_cent.get(Symmetry{d})), ss, coor,
+                                                    TensorDimens(sym, true_nvs));
+              this->insert(std::move(ten));
             }
         }
     }
@@ -439,9 +441,9 @@ public:
   }
 protected:
   template <int t>
-  _Ttensor *faaDiBrunoZ(const Symmetry &sym) const;
+  std::unique_ptr<_Ttensor> faaDiBrunoZ(const Symmetry &sym) const;
   template <int t>
-  _Ttensor *faaDiBrunoG(const Symmetry &sym) const;
+  std::unique_ptr<_Ttensor> faaDiBrunoG(const Symmetry &sym) const;
 
   // convenience access methods
   template<int t>
@@ -472,12 +474,12 @@ protected:
    of a given symmetry. */
 
 template <int t>
-_Ttensor *
+std::unique_ptr<_Ttensor>
 KOrderStoch::faaDiBrunoZ(const Symmetry &sym) const
 {
   JournalRecordPair pa(journal);
   pa << "Faa Di Bruno ZX container for " << sym << endrec;
-  _Ttensor *res = new _Ttensor(ypart.ny(), TensorDimens(sym, nvs));
+  auto res = std::make_unique<_Ttensor>(ypart.ny(), TensorDimens(sym, nvs));
   FaaDiBruno bruno(journal);
   bruno.calculate(Zstack<t>(), f, *res);
   return res;
@@ -487,13 +489,13 @@ KOrderStoch::faaDiBrunoZ(const Symmetry &sym) const
    $G(y,u,\sigma)=h(g^*(y,u,\sigma),\sigma)$ of a given symmetry. */
 
 template <int t>
-_Ttensor *
+std::unique_ptr<_Ttensor>
 KOrderStoch::faaDiBrunoG(const Symmetry &sym) const
 {
   JournalRecordPair pa(journal);
   pa << "Faa Di Bruno GX container for " << sym << endrec;
   TensorDimens tdims(sym, nvs);
-  auto *res = new _Ttensor(ypart.nyss(), tdims);
+  auto res = std::make_unique<_Ttensor>(ypart.nyss(), tdims);
   FaaDiBruno bruno(journal);
   bruno.calculate(Gstack<t>(), h<t>(), *res);
   return res;
@@ -524,16 +526,18 @@ KOrderStoch::performStep(int order)
           JournalRecordPair pa(journal);
           pa << "Recovering symmetry " << si << endrec;
 
-          _Ttensor *G_sym = faaDiBrunoG<t>(si);
-          G<t>().insert(G_sym);
+          auto G_sym = faaDiBrunoG<t>(si);
+          auto G_sym_ptr = G_sym.get();
+          G<t>().insert(std::move(G_sym));
 
-          _Ttensor *g_sym = faaDiBrunoZ<t>(si);
+          auto g_sym = faaDiBrunoZ<t>(si);
+          auto g_sym_ptr = g_sym.get();
           g_sym->mult(-1.0);
           matA.multInv(*g_sym);
-          g<t>().insert(g_sym);
-          gs<t>().insert(new _Ttensor(ypart.nstat, ypart.nys(), *g_sym));
+          g<t>().insert(std::move(g_sym));
+          gs<t>().insert(std::make_unique<_Ttensor>(ypart.nstat, ypart.nys(), *g_sym_ptr));
 
-          Gstack<t>().multAndAdd(1, h<t>(), *G_sym);
+          Gstack<t>().multAndAdd(1, h<t>(), *G_sym_ptr);
         }
     }
 }
diff --git a/dynare++/kord/tests.cc b/dynare++/kord/tests.cc
index dbdf72c23..10e3483e1 100644
--- a/dynare++/kord/tests.cc
+++ b/dynare++/kord/tests.cc
@@ -63,17 +63,17 @@ Rand::discrete(double prob)
 
 struct SparseGenerator
 {
-  static FSSparseTensor *makeTensor(int dim, int nv, int r,
-                                    double fill, double m);
+  static std::unique_ptr<FSSparseTensor> makeTensor(int dim, int nv, int r,
+                                                    double fill, double m);
   static void fillContainer(TensorContainer<FSSparseTensor> &c,
                             int maxdim, int nv, int r, double m);
 };
 
-FSSparseTensor *
+std::unique_ptr<FSSparseTensor>
 SparseGenerator::makeTensor(int dim, int nv, int r,
                             double fill, double m)
 {
-  auto *res = new FSSparseTensor(dim, nv, r);
+  auto res = std::make_unique<FSSparseTensor>(dim, nv, r);
   FFSTensor dummy(0, nv, dim);
   for (Tensor::index fi = dummy.begin(); fi != dummy.end(); ++fi)
     for (int i = 0; i < r; i++)
diff --git a/dynare++/src/dynare3.cc b/dynare++/src/dynare3.cc
index 279c31003..2b5f85942 100644
--- a/dynare++/src/dynare3.cc
+++ b/dynare++/src/dynare3.cc
@@ -345,8 +345,8 @@ DynareDerEvalLoader::DynareDerEvalLoader(const ogp::FineAtoms &a,
   md.clear();
   for (int iord = 1; iord <= order; iord++)
     {
-      auto *t = new FSSparseTensor(iord, atoms.ny()+atoms.nys()+atoms.nyss()+atoms.nexo(), atoms.ny());
-      md.insert(t);
+      auto t = std::make_unique<FSSparseTensor>(iord, atoms.ny()+atoms.nys()+atoms.nyss()+atoms.nexo(), atoms.ny());
+      md.insert(std::move(t));
     }
 }
 
diff --git a/dynare++/tl/cc/normal_moments.cc b/dynare++/tl/cc/normal_moments.cc
index 3f6f96853..f0af6968c 100644
--- a/dynare++/tl/cc/normal_moments.cc
+++ b/dynare++/tl/cc/normal_moments.cc
@@ -29,9 +29,9 @@ UNormalMoments::generateMoments(int maxdim, const TwoDMatrix &v)
               "Variance-covariance matrix is not square in UNormalMoments constructor");
 
   int nv = v.nrows();
-  auto *mom2 = new URSingleTensor(nv, 2);
+  auto mom2 = std::make_unique<URSingleTensor>(nv, 2);
   mom2->getData() = v.getData();
-  insert(mom2);
+  insert(std::move(mom2));
   auto kronv = std::make_unique<URSingleTensor>(nv, 2);
   kronv->getData() = v.getData();
   for (int d = 4; d <= maxdim; d += 2)
@@ -41,7 +41,7 @@ UNormalMoments::generateMoments(int maxdim, const TwoDMatrix &v)
                          ConstVector(kronv->getData()),
                          newkronv->getData());
       kronv = std::move(newkronv);
-      auto *mom = new URSingleTensor(nv, d);
+      auto mom = std::make_unique<URSingleTensor>(nv, d);
       // apply $F_n$ to |kronv|
       /* Here we go through all equivalences, select only those having 2
          elements in each class, then go through all elements in |kronv| and
@@ -66,7 +66,7 @@ UNormalMoments::generateMoments(int maxdim, const TwoDMatrix &v)
                 mom->get(*it2, 0) += kronv->get(*it, 0);
               }
           }
-      insert(mom);
+      insert(std::move(mom));
     }
 }
 
@@ -89,8 +89,5 @@ FNormalMoments::FNormalMoments(const UNormalMoments &moms)
   : TensorContainer<FRSingleTensor>(1)
 {
   for (const auto &mom : moms)
-    {
-      auto *fm = new FRSingleTensor(*(mom.second));
-      insert(fm);
-    }
+    insert(std::make_unique<FRSingleTensor>(*(mom.second)));
 }
diff --git a/dynare++/tl/cc/t_container.cc b/dynare++/tl/cc/t_container.cc
index d9202d66d..720c25ca3 100644
--- a/dynare++/tl/cc/t_container.cc
+++ b/dynare++/tl/cc/t_container.cc
@@ -11,11 +11,8 @@ const int FGSContainer::num_one_time = 10;
 UGSContainer::UGSContainer(const FGSContainer &c)
   : TensorContainer<UGSTensor>(c.num())
 {
-  for (const auto & it : c)
-    {
-      auto *unfolded = new UGSTensor(*(it.second));
-      insert(unfolded);
-    }
+  for (const auto &it : c)
+    insert(std::make_unique<UGSTensor>(*(it.second)));
 }
 
 /* We set |l| to dimension of |t|, this is a tensor which multiplies
@@ -58,10 +55,7 @@ FGSContainer::FGSContainer(const UGSContainer &c)
   : TensorContainer<FGSTensor>(c.num())
 {
   for (const auto & it : c)
-    {
-      auto *folded = new FGSTensor(*(it.second));
-      insert(folded);
-    }
+    insert(std::make_unique<FGSTensor>(*(it.second)));
 }
 
 // |FGSContainer::multAndAdd| folded code
diff --git a/dynare++/tl/cc/t_container.hh b/dynare++/tl/cc/t_container.hh
index a467e5ad4..7f39749fb 100644
--- a/dynare++/tl/cc/t_container.hh
+++ b/dynare++/tl/cc/t_container.hh
@@ -59,6 +59,8 @@
 #include <map>
 #include <string>
 #include <sstream>
+#include <memory>
+#include <utility>
 
 #include <matio.h>
 
@@ -96,8 +98,7 @@ class TensorContainer
 protected:
   using _const_ptr = const _Ttype *;
   using _ptr = _Ttype *;
-  using _Map = std::map<Symmetry, _ptr, ltsym>;
-  using _mvtype = typename _Map::value_type;
+  using _Map = std::map<Symmetry, std::unique_ptr<_Ttype>, ltsym>;
 public:
   using iterator = typename _Map::iterator;
   using const_iterator = typename _Map::const_iterator;
@@ -113,14 +114,12 @@ public:
   }
   /* This is just a copy constructor. This makes a hard copy of all tensors. */
   TensorContainer(const TensorContainer<_Ttype> &c)
-    : n(c.n), m(), ebundle(c.ebundle)
+    : n(c.n), ebundle(c.ebundle)
   {
-    for (auto it = c.m.begin(); it != c.m.end(); ++it)
-      {
-        auto *ten = new _Ttype(*((*it).second));
-        insert(ten);
-      }
+    for (const auto &it : c.m)
+      insert(std::make_unique<_Ttype>(*(it.second)));
   }
+  TensorContainer(TensorContainer<_Ttype> &&) = default;
 
   // |TensorContainer| subtensor constructor
   /* This constructor constructs a new tensor container, whose tensors
@@ -128,11 +127,8 @@ public:
   TensorContainer(int first_row, int num, TensorContainer<_Ttype> &c)
     : n(c.n), ebundle(*(tls.ebundle))
   {
-    for (auto it = c.m.begin(); it != c.m.end(); ++it)
-      {
-        auto *t = new _Ttype(first_row, num, *((*it).second));
-        insert(t);
-      }
+    for (const auto &it : c.m)
+      insert(std::make_unique<_Ttype>(first_row, num, *(it.second)));
   }
 
   _const_ptr
@@ -147,9 +143,7 @@ public:
         return nullptr;
       }
     else
-      {
-        return (*it).second;
-      }
+      return it->second.get();
   }
 
   _ptr
@@ -164,9 +158,7 @@ public:
         return nullptr;
       }
     else
-      {
-        return (*it).second;
-      }
+      return it->second.get();
   }
 
   bool
@@ -178,49 +170,37 @@ public:
     return it != m.end();
   }
 
-  void
-  insert(_ptr t)
+  virtual void
+  insert(std::unique_ptr<_Ttype> t)
   {
     TL_RAISE_IF(t->getSym().num() != num(),
                 "Incompatible symmetry insertion in TensorContainer::insert");
     TL_RAISE_IF(check(t->getSym()),
                 "Tensor already in container in TensorContainer::insert");
-    m.insert(_mvtype(t->getSym(), t));
     if (!t->isFinite())
-      {
-        throw TLException(__FILE__, __LINE__,  "NaN or Inf asserted in TensorContainer::insert");
-      }
+      throw TLException(__FILE__, __LINE__,  "NaN or Inf asserted in TensorContainer::insert");
+    m.emplace(t->getSym(), std::move(t));
   }
 
   void
   remove(const Symmetry &s)
   {
-    auto it = m.find(s);
-    if (it != m.end())
-      {
-        _ptr t = (*it).second;
-        m.erase(it);
-        delete t;
-      }
+    m.erase(s);
   }
 
   void
   clear()
   {
-    while (!m.empty())
-      {
-        delete (*(m.begin())).second;
-        m.erase(m.begin());
-      }
+    m.clear();
   }
 
   int
   getMaxDim() const
   {
     int res = -1;
-    for (auto run = m.begin(); run != m.end(); ++run)
+    for (const auto &run : m)
       {
-        int dim = (*run).first.dimen();
+        int dim = run.first.dimen();
         if (dim > res)
           res = dim;
       }
@@ -235,8 +215,8 @@ public:
     for (const_iterator it = m.begin(); it != m.end(); ++it)
       {
         printf("Symmetry: ");
-        (*it).first.print();
-        ((*it).second)->print();
+        (it->first).print();
+        (it->second)->print();
       }
   }
 
@@ -255,7 +235,7 @@ public:
             sprintf(tmp, "_%d", sym[i]);
             strcat(lname, tmp);
           }
-        ConstTwoDMatrix m(*((*it).second));
+        ConstTwoDMatrix m(*(it->second));
         m.writeMat(fd, lname);
       }
   }
@@ -294,10 +274,7 @@ public:
     return res;
   }
 
-  virtual ~TensorContainer()
-  {
-    clear();
-  }
+  virtual ~TensorContainer() = default;
 
   int
   num() const
@@ -342,9 +319,8 @@ public:
     : TensorContainer<UGSTensor>(nn)
   {
   }
-  UGSContainer(const UGSContainer &uc)
-     
-  = default;
+  UGSContainer(const UGSContainer &) = default;
+  UGSContainer(UGSContainer &&) = default;
   UGSContainer(const FGSContainer &c);
   void multAndAdd(const UGSTensor &t, UGSTensor &out) const;
 };
@@ -370,9 +346,8 @@ public:
     : TensorContainer<FGSTensor>(nn)
   {
   }
-  FGSContainer(const FGSContainer &fc)
-     
-  = default;
+  FGSContainer(const FGSContainer &) = default;
+  FGSContainer(FGSContainer &&) = default;
   FGSContainer(const UGSContainer &c);
   void multAndAdd(const FGSTensor &t, FGSTensor &out) const;
   void multAndAdd(const UGSTensor &t, FGSTensor &out) const;
diff --git a/dynare++/tl/cc/t_polynomial.cc b/dynare++/tl/cc/t_polynomial.cc
index 941e29584..f1155d211 100644
--- a/dynare++/tl/cc/t_polynomial.cc
+++ b/dynare++/tl/cc/t_polynomial.cc
@@ -50,17 +50,13 @@ PowerProvider::~PowerProvider()
 UTensorPolynomial::UTensorPolynomial(const FTensorPolynomial &fp)
   : TensorPolynomial<UFSTensor, UGSTensor, URSingleTensor>(fp.nrows(), fp.nvars())
 {
-  for (const auto & it : fp)
-    {
-      insert(new UFSTensor(*(it.second)));
-    }
+  for (const auto &it : fp)
+    insert(std::make_unique<UFSTensor>(*(it.second)));
 }
 
 FTensorPolynomial::FTensorPolynomial(const UTensorPolynomial &up)
   : TensorPolynomial<FFSTensor, FGSTensor, FRSingleTensor>(up.nrows(), up.nvars())
 {
-  for (const auto & it : up)
-    {
-      insert(new FFSTensor(*(it.second)));
-    }
+  for (const auto &it : up)
+    insert(std::make_unique<FFSTensor>(*(it.second)));
 }
diff --git a/dynare++/tl/cc/t_polynomial.hh b/dynare++/tl/cc/t_polynomial.hh
index cefd105ba..84a259b48 100644
--- a/dynare++/tl/cc/t_polynomial.hh
+++ b/dynare++/tl/cc/t_polynomial.hh
@@ -175,14 +175,13 @@ public:
                 /* The pointer |ten| is either a new tensor or got from |this| container. */
                 _Ttype *ten;
                 if (_Tparent::check(Symmetry{j}))
-                  {
-                    ten = _Tparent::get(Symmetry{j});
-                  }
+                  ten = _Tparent::get(Symmetry{j});
                 else
                   {
-                    ten = new _Ttype(nrows(), nvars(), j);
-                    ten->zeros();
-                    insert(ten);
+                    auto ten_smart = std::make_unique<_Ttype>(nrows(), nvars(), j);
+                    ten_smart->zeros();
+                    ten = ten_smart.get();
+                    insert(std::move(ten_smart));
                   }
 
                 Symmetry sym{i, j};
@@ -208,14 +207,13 @@ public:
             /* Same code as above */
             _Ttype *ten;
             if (_Tparent::check(Symmetry{j}))
-              {
-                ten = _Tparent::get(Symmetry{j});
-              }
+              ten = _Tparent::get(Symmetry{j});
             else
               {
-                ten = new _Ttype(nrows(), nvars(), j);
-                ten->zeros();
-                insert(ten);
+                auto ten_smart = std::make_unique<_Ttype>(nrows(), nvars(), j);
+                ten_smart->zeros();
+                ten = ten_smart.get();
+                insert(std::move(ten_smart));
               }
 
             Symmetry sym{0, j};
@@ -308,15 +306,15 @@ public:
      number of variables. Then we insert and update the |maxdim|. */
 
   void
-  insert(_ptr t)
+  insert(std::unique_ptr<_Ttype> t) override
   {
     TL_RAISE_IF(t->nrows() != nr,
                 "Wrong number of rows in TensorPolynomial::insert");
     TL_RAISE_IF(t->nvar() != nv,
                 "Wrong number of variables in TensorPolynomial::insert");
-    TensorContainer<_Ttype>::insert(t);
     if (maxdim < t->dimen())
       maxdim = t->dimen();
+    TensorContainer<_Ttype>::insert(std::move(t));
   }
 
   /* The polynomial takes the form
@@ -359,13 +357,13 @@ public:
      |g.derivative(0)|, |g.derivative(1)| and |der=g.evalPartially(2, v)|
      calculates $2!$ multiple of the second derivative of |g| at |v|. */
 
-  _Ttype *
+  std::unique_ptr<_Ttype>
   evalPartially(int s, const ConstVector &v)
   {
     TL_RAISE_IF(v.length() != nvars(),
                 "Wrong length of vector for TensorPolynomial::evalPartially");
 
-    auto *res = new _Ttype(nrows(), nvars(), s);
+    auto res = std::make_unique<_Ttype>(nrows(), nvars(), s);
     res->zeros();
 
     if (_Tparent::check(Symmetry{s}))
@@ -376,15 +374,13 @@ public:
         if (_Tparent::check(Symmetry{d}))
           {
             const _Ttype &ltmp = *(_Tparent::get(Symmetry{d}));
-            auto *last = new _Ttype(ltmp);
+            auto last = std::make_unique<_Ttype>(ltmp);
             for (int j = 0; j < d - s; j++)
               {
-                auto *newlast = new _Ttype(*last, v);
-                delete last;
-                last = newlast;
+                auto newlast = std::make_unique<_Ttype>(*last, v);
+                last = std::move(newlast);
               }
             res->add(1.0, *last);
-            delete last;
           }
       }
 
diff --git a/dynare++/tl/testing/factory.hh b/dynare++/tl/testing/factory.hh
index dca4e5009..284dc6f28 100644
--- a/dynare++/tl/testing/factory.hh
+++ b/dynare++/tl/testing/factory.hh
@@ -5,6 +5,7 @@
 #define FACTORY_H
 
 #include <random>
+#include <memory>
 
 #include "symmetry.hh"
 #include "int_sequence.hh"
@@ -25,10 +26,10 @@ public:
   double get();
   // this can be used with UGSTensor, FGSTensor
   template <class _Ttype>
-  _Ttype *
+  std::unique_ptr<_Ttype>
   make(int r, const Symmetry &s, const IntSequence &nvs)
   {
-    _Ttype *res = new _Ttype(r, TensorDimens(s, nvs));
+    auto res = std::make_unique<_Ttype>(r, TensorDimens(s, nvs));
     init(s, nvs);
     fillMatrix(*res);
     return res;
@@ -36,10 +37,10 @@ public:
 
   // this can be used with FFSTensor, UFSTensor, FRTensor, URTensor
   template <class _Ttype>
-  _Ttype *
+  std::unique_ptr<_Ttype>
   make(int r, int nv, int dim)
   {
-    auto *res = new _Ttype(r, nv, dim);
+    auto res = std::make_unique<_Ttype>(r, nv, dim);
     init(dim, nv);
     fillMatrix(*res);
     return res;
@@ -56,8 +57,7 @@ public:
         {
           // full symmetry
           Symmetry sym{dim};
-          auto *t = make<_Ttype>(r, sym, nvs);
-          res->insert(t);
+          res->insert(make<_Ttype>(r, sym, nvs));
         }
       else
         {
@@ -65,8 +65,7 @@ public:
           for (int i = 0; i <= dim; i++)
             {
               Symmetry sym{i, dim-i};
-              auto *t = make<_Ttype>(r, sym, nvs);
-              res->insert(t);
+              res->insert(make<_Ttype>(r, sym, nvs));
             }
         }
     return res;
@@ -78,10 +77,7 @@ public:
   {
     auto *p = new _Ptype(r, nv);
     for (int d = 1; d <= maxdim; d++)
-      {
-        auto *t = make<_Ttype>(r, nv, d);
-        p->insert(t);
-      }
+      p->insert(make<_Ttype>(r, nv, d));
     return p;
   }
 
diff --git a/dynare++/tl/testing/monoms.cc b/dynare++/tl/testing/monoms.cc
index 737fed31c..89c33308f 100644
--- a/dynare++/tl/testing/monoms.cc
+++ b/dynare++/tl/testing/monoms.cc
@@ -102,11 +102,11 @@ Monom1Vector::deriv(const IntSequence &c, Vector &out) const
     out[i] = x[i].deriv(c);
 }
 
-FGSTensor *
+std::unique_ptr<FGSTensor>
 Monom1Vector::deriv(int dim) const
 {
-  FGSTensor *res
-    = new FGSTensor(len, TensorDimens(Symmetry{dim}, IntSequence(1, nx)));
+  auto res = std::make_unique<FGSTensor>(len, TensorDimens(Symmetry{dim},
+                                                           IntSequence(1, nx)));
   for (Tensor::index it = res->begin(); it != res->end(); ++it)
     {
       Vector outcol{res->getCol(*it)};
@@ -176,11 +176,11 @@ Monom2Vector::deriv(const Symmetry &s, const IntSequence &c,
     out[i] = y[i].deriv(cy) * u[i].deriv(cu);
 }
 
-FGSTensor *
+std::unique_ptr<FGSTensor>
 Monom2Vector::deriv(const Symmetry &s) const
 {
   IntSequence nvs{ny, nu};
-  FGSTensor *t = new FGSTensor(len, TensorDimens(s, nvs));
+  auto t = std::make_unique<FGSTensor>(len, TensorDimens(s, nvs));
   for (Tensor::index it = t->begin(); it != t->end(); ++it)
     {
       Vector col{t->getCol(*it)};
@@ -318,11 +318,11 @@ Monom4Vector::deriv(const Symmetry &s, const IntSequence &coor,
     }
 }
 
-FGSTensor *
+std::unique_ptr<FGSTensor>
 Monom4Vector::deriv(const Symmetry &s) const
 {
   IntSequence nvs{nx1, nx2, nx3, nx4};
-  FGSTensor *res = new FGSTensor(len, TensorDimens(s, nvs));
+  auto res = std::make_unique<FGSTensor>(len, TensorDimens(s, nvs));
   for (Tensor::index run = res->begin(); run != res->end(); ++run)
     {
       Vector col{res->getCol(*run)};
@@ -416,8 +416,7 @@ SparseDerivGenerator::~SparseDerivGenerator()
 
 DenseDerivGenerator::DenseDerivGenerator(int ng, int nx, int ny, int nu,
                                          int mx, double prob, int maxdim)
-  : maxdimen(maxdim), ts(new FGSTensor *[maxdimen]),
-    uts(new UGSTensor *[maxdimen])
+  : maxdimen(maxdim), ts(maxdimen), uts(maxdimen)
 {
   intgen.init(ng, nx, ny, nu, nu, mx, prob);
   Monom1Vector g(nx, ng);
@@ -427,10 +426,7 @@ DenseDerivGenerator::DenseDerivGenerator(int ng, int nx, int ny, int nu,
   rcont = r.deriv(maxdimen);
   uxcont = nullptr;
   for (int d = 1; d <= maxdimen; d++)
-    {
-      ts[d-1] = g.deriv(d);
-      uts[d-1] = nullptr;
-    }
+    ts[d-1] = g.deriv(d);
 }
 
 void
@@ -438,19 +434,11 @@ DenseDerivGenerator::unfold()
 {
   uxcont = new UGSContainer(*xcont);
   for (int i = 0; i < maxdimen; i++)
-    uts[i] = new UGSTensor(*(ts[i]));
+    uts[i] = std::make_unique<UGSTensor>(*(ts[i]));
 }
 
 DenseDerivGenerator::~DenseDerivGenerator()
 {
   delete xcont;
   delete rcont;
-  for (int i = 0; i < maxdimen; i++)
-    {
-      delete ts[i];
-      if (uts[i])
-        delete uts[i];
-    }
-  delete [] ts;
-  delete [] uts;
 }
diff --git a/dynare++/tl/testing/monoms.hh b/dynare++/tl/testing/monoms.hh
index d6716bfd5..b38e4c925 100644
--- a/dynare++/tl/testing/monoms.hh
+++ b/dynare++/tl/testing/monoms.hh
@@ -6,6 +6,7 @@
 
 #include <vector>
 #include <random>
+#include <memory>
 
 #include "int_sequence.hh"
 #include "gs_tensor.hh"
@@ -49,7 +50,7 @@ public:
   Monom1Vector(int nxx, int l);
   ~Monom1Vector() = default;
   void deriv(const IntSequence &c, Vector &out) const;
-  FGSTensor *deriv(int dim) const;
+  std::unique_ptr<FGSTensor> deriv(int dim) const;
   void print() const;
 };
 
@@ -66,7 +67,7 @@ public:
   Monom2Vector(const Monom1Vector &g, const Monom2Vector &xmon);
   ~Monom2Vector() = default;
   void deriv(const Symmetry &s, const IntSequence &c, Vector &out) const;
-  FGSTensor *deriv(const Symmetry &s) const;
+  std::unique_ptr<FGSTensor> deriv(const Symmetry &s) const;
   FGSContainer *deriv(int maxdim) const;
   void print() const;
 };
@@ -88,7 +89,7 @@ public:
                const Monom4Vector &g);
   ~Monom4Vector() = default;
   FSSparseTensor *deriv(int dim) const;
-  FGSTensor *deriv(const Symmetry &s) const;
+  std::unique_ptr<FGSTensor> deriv(const Symmetry &s) const;
   void deriv(const Symmetry &s, const IntSequence &coor, Vector &out) const;
   void print() const;
 protected:
@@ -112,9 +113,9 @@ struct DenseDerivGenerator
   int maxdimen;
   FGSContainer *xcont;
   FGSContainer *rcont;
-  FGSTensor **const ts;
+  std::vector<std::unique_ptr<FGSTensor>> ts;
   UGSContainer *uxcont;
-  UGSTensor **const uts;
+  std::vector<std::unique_ptr<UGSTensor>> uts;
   DenseDerivGenerator(int ng, int nx, int ny, int nu,
                       int mx, double prob, int maxdim);
   void unfold();
diff --git a/dynare++/tl/testing/tests.cc b/dynare++/tl/testing/tests.cc
index 47635a597..6bbb97f85 100644
--- a/dynare++/tl/testing/tests.cc
+++ b/dynare++/tl/testing/tests.cc
@@ -45,27 +45,24 @@ protected:
   template <class _Ttype>
   static bool index_offset(const Symmetry &s, const IntSequence &nvs);
 
-  static bool fold_unfold(const FTensor *folded);
+  static bool fold_unfold(std::unique_ptr<FTensor> folded);
   static bool
   fs_fold_unfold(int r, int nv, int dim)
   {
     Factory f;
-    FTensor *folded = f.make<FFSTensor>(r, nv, dim);
-    return fold_unfold(folded); // folded deallocated in fold_unfold
+    return fold_unfold(f.make<FFSTensor>(r, nv, dim));
   }
   static bool
   r_fold_unfold(int r, int nv, int dim)
   {
     Factory f;
-    FTensor *folded = f.make<FRTensor>(r, nv, dim);
-    return fold_unfold(folded); // folded deallocated in fold_unfold
+    return fold_unfold(f.make<FRTensor>(r, nv, dim));
   }
   static bool
   gs_fold_unfold(int r, const Symmetry &s, const IntSequence &nvs)
   {
     Factory f;
-    FTensor *folded = f.make<FGSTensor>(r, s, nvs);
-    return fold_unfold(folded); // folded deallocated in fold_unfold
+    return fold_unfold(f.make<FGSTensor>(r, s, nvs));
   }
 
   static bool dense_prod(const Symmetry &bsym, const IntSequence &bnvs,
@@ -194,7 +191,7 @@ TestRunnable::index_offset(const Symmetry &s, const IntSequence &nvs)
 }
 
 bool
-TestRunnable::fold_unfold(const FTensor *folded)
+TestRunnable::fold_unfold(std::unique_ptr<FTensor> folded)
 {
   auto unfolded = folded->unfold();
   auto folded2 = unfolded->fold();
@@ -206,8 +203,6 @@ TestRunnable::fold_unfold(const FTensor *folded)
             << "\tdifference normInf: " << normInf << '\n'
             << "\tdifference norm1:   " << norm1 << '\n';
 
-  delete folded;
-
   return normInf < 1.0e-15;
 }
 
@@ -218,8 +213,7 @@ TestRunnable::dense_prod(const Symmetry &bsym, const IntSequence &bnvs,
   Factory f;
   FGSContainer *cont
     = f.makeCont<FGSTensor, FGSContainer>(hnv, bnvs, bsym.dimen()-hdim+1);
-  auto *fh
-    = f.make<FGSTensor>(rows, Symmetry{hdim}, IntSequence(1, hnv));
+  auto fh = f.make<FGSTensor>(rows, Symmetry{hdim}, IntSequence(1, hnv));
   UGSTensor uh(*fh);
   FGSTensor fb(rows, TensorDimens(bsym, bnvs));
   fb.getData().zeros();
@@ -248,7 +242,6 @@ TestRunnable::dense_prod(const Symmetry &bsym, const IntSequence &bnvs,
             << "\tunfolded difference normInf: " << normInf << '\n';
 
   delete cont;
-  delete fh;
 
   return norm < 1.e-13;
 }
@@ -431,7 +424,7 @@ TestRunnable::folded_contraction(int r, int nv, int dim)
   Factory fact;
   Vector x{fact.makeVector(nv)};
 
-  auto *forig = fact.make<FFSTensor>(r, nv, dim);
+  auto forig = fact.make<FFSTensor>(r, nv, dim);
   auto *f = new FFSTensor(*forig);
   clock_t ctime = clock();
   for (int d = dim-1; d > 0; d--)
@@ -472,9 +465,8 @@ TestRunnable::unfolded_contraction(int r, int nv, int dim)
   Factory fact;
   Vector x{fact.makeVector(nv)};
 
-  auto *forig = fact.make<FFSTensor>(r, nv, dim);
+  auto forig = fact.make<FFSTensor>(r, nv, dim);
   UFSTensor uorig(*forig);
-  delete forig;
   auto *u = new UFSTensor(uorig);
   clock_t ctime = clock();
   for (int d = dim-1; d > 0; d--)
diff --git a/mex/sources/k_order_perturbation/k_ord_dynare.cc b/mex/sources/k_order_perturbation/k_ord_dynare.cc
index 12827e5fe..700f61fc6 100644
--- a/mex/sources/k_order_perturbation/k_ord_dynare.cc
+++ b/mex/sources/k_order_perturbation/k_ord_dynare.cc
@@ -51,7 +51,7 @@ KordpDynare::KordpDynare(const std::vector<std::string> &endo, int num_endo,
 
   //	Initialise ModelDerivativeContainer(*this, this->md, nOrder);
   for (int iord = 1; iord <= nOrder; iord++)
-    md.insert(new FSSparseTensor(iord, nY+nYs+nYss+nExog, nY));
+    md.insert(std::make_unique<FSSparseTensor>(iord, nY+nYs+nYss+nExog, nY));
 }
 
 KordpDynare::KordpDynare(const std::vector<std::string> &endo, int num_endo,
@@ -73,7 +73,7 @@ KordpDynare::KordpDynare(const std::vector<std::string> &endo, int num_endo,
 
   //	Initialise ModelDerivativeContainer(*this, this->md, nOrder);
   for (int iord = 1; iord <= nOrder; iord++)
-    md.insert(new FSSparseTensor(iord, nY+nYs+nYss+nExog, nY));
+    md.insert(std::make_unique<FSSparseTensor>(iord, nY+nYs+nYss+nExog, nY));
 }
 
 KordpDynare::~KordpDynare()
@@ -162,7 +162,7 @@ void
 KordpDynare::populateDerivativesContainer(const TwoDMatrix &g, int ord, const std::vector<int> &vOrder)
 {
   // model derivatives FSSparseTensor instance
-  FSSparseTensor *mdTi = (new FSSparseTensor(ord, nJcols, nY));
+  auto mdTi = std::make_unique<FSSparseTensor>(ord, nJcols, nY);
 
   IntSequence s(ord, 0);
 
@@ -249,8 +249,7 @@ KordpDynare::populateDerivativesContainer(const TwoDMatrix &g, int ord, const st
 
   // md container
   md.remove(Symmetry{ord});
-  md.insert(mdTi);
-  // No need to delete mdTi, it will be deleted by TensorContainer destructor
+  md.insert(std::move(mdTi));
 }
 
 /*********************************************************