diff --git a/dynare++/extern/matlab/dynare_simul.cc b/dynare++/extern/matlab/dynare_simul.cc
index 56486d608..b36d4f10a 100644
--- a/dynare++/extern/matlab/dynare_simul.cc
+++ b/dynare++/extern/matlab/dynare_simul.cc
@@ -26,6 +26,8 @@
 #include "fs_tensor.hh"
 #include "SylvException.hh"
 
+#include <string>
+
 extern "C" {
   void
   mexFunction(int nlhs, mxArray *plhs[],
@@ -34,20 +36,20 @@ extern "C" {
     if (nhrs < 12 || nlhs != 2)
       DYN_MEX_FUNC_ERR_MSG_TXT("dynare_simul_ must have at least 12 input parameters and exactly 2 output arguments.\n");
 
-    auto order = (int) mxGetScalar(prhs[0]);
+    int order = static_cast<int>(mxGetScalar(prhs[0]));
     if (nhrs != 12 + order)
       DYN_MEX_FUNC_ERR_MSG_TXT("dynare_simul_ must have exactly 11+order input parameters.\n");
 
-    auto nstat = (int) mxGetScalar(prhs[1]);
-    auto npred = (int) mxGetScalar(prhs[2]);
-    auto nboth = (int) mxGetScalar(prhs[3]);
-    auto nforw = (int) mxGetScalar(prhs[4]);
-    auto nexog = (int) mxGetScalar(prhs[5]);
+    int nstat = static_cast<int>(mxGetScalar(prhs[1]));
+    int npred = static_cast<int>(mxGetScalar(prhs[2]));
+    int nboth = static_cast<int>(mxGetScalar(prhs[3]));
+    int nforw = static_cast<int>(mxGetScalar(prhs[4]));
+    int nexog = static_cast<int>(mxGetScalar(prhs[5]));
 
     const mxArray *const ystart = prhs[6];
     const mxArray *const shocks = prhs[7];
     const mxArray *const vcov = prhs[8];
-    auto seed = (int) mxGetScalar(prhs[9]);
+    int seed = static_cast<int>(mxGetScalar(prhs[9]));
     const mxArray *const ysteady = prhs[10];
     const mwSize *const ystart_dim = mxGetDimensions(ystart);
     const mwSize *const shocks_dim = mxGetDimensions(shocks);
@@ -55,18 +57,18 @@ extern "C" {
     const mwSize *const ysteady_dim = mxGetDimensions(ysteady);
 
     int ny = nstat + npred + nboth + nforw;
-    if (ny != (int) ystart_dim[0])
+    if (ny != static_cast<int>(ystart_dim[0]))
       DYN_MEX_FUNC_ERR_MSG_TXT("ystart has wrong number of rows.\n");
     if (1 != ystart_dim[1])
       DYN_MEX_FUNC_ERR_MSG_TXT("ystart has wrong number of cols.\n");
     int nper = shocks_dim[1];
-    if (nexog != (int) shocks_dim[0])
+    if (nexog != static_cast<int>(shocks_dim[0]))
       DYN_MEX_FUNC_ERR_MSG_TXT("shocks has a wrong number of rows.\n");
-    if (nexog != (int) vcov_dim[0])
+    if (nexog != static_cast<int>(vcov_dim[0]))
       DYN_MEX_FUNC_ERR_MSG_TXT("vcov has a wrong number of rows.\n");
-    if (nexog != (int) vcov_dim[1])
+    if (nexog != static_cast<int>(vcov_dim[1]))
       DYN_MEX_FUNC_ERR_MSG_TXT("vcov has a wrong number of cols.\n");
-    if (ny != (int) ysteady_dim[0])
+    if (ny != static_cast<int>(ysteady_dim[0]))
       DYN_MEX_FUNC_ERR_MSG_TXT("ysteady has wrong number of rows.\n");
     if (1 != ysteady_dim[1])
       DYN_MEX_FUNC_ERR_MSG_TXT("ysteady has wrong number of cols.\n");
@@ -85,29 +87,18 @@ extern "C" {
             const mxArray *gk = prhs[11+dim];
             const mwSize *const gk_dim = mxGetDimensions(gk);
             FFSTensor ft(ny, npred+nboth+nexog, dim);
-            if (ft.ncols() != (int) gk_dim[1])
-              {
-                char buf[1000];
-                sprintf(buf, "Wrong number of columns for folded tensor: got %d but I want %d\n",
-                        (int) gk_dim[1], ft.ncols());
-                DYN_MEX_FUNC_ERR_MSG_TXT(buf);
-              }
-            if (ft.nrows() != (int) gk_dim[0])
-              {
-                char buf[1000];
-                sprintf(buf, "Wrong number of rows for folded tensor: got %d but I want %d\n",
-                        (int) gk_dim[0], ft.nrows());
-                DYN_MEX_FUNC_ERR_MSG_TXT(buf);
-              }
+            if (ft.ncols() != static_cast<int>(gk_dim[1]))
+              DYN_MEX_FUNC_ERR_MSG_TXT("Wrong number of columns for folded tensor: got " + std::to_string(gk_dim[1]) + " but i want " + std::to_string(ft.ncols()) + '\n');
+            if (ft.nrows() != static_cast<int>(gk_dim[0]))
+              DYN_MEX_FUNC_ERR_MSG_TXT("Wrong number of rows for folded tensor: got " + std::to_string(gk_dim[0]) + " but i want " + std::to_string(ft.nrows()) + '\n');
             ft.zeros();
             ConstTwoDMatrix gk_mat(ft.nrows(), ft.ncols(), ConstVector{gk});
             ft.add(1.0, gk_mat);
             pol.insert(std::make_unique<UFSTensor>(ft));
           }
         // form the decision rule
-        UnfoldDecisionRule
-          dr(pol, PartitionY(nstat, npred, nboth, nforw),
-             nexog, ConstVector{ysteady});
+        UnfoldDecisionRule dr(pol, PartitionY(nstat, npred, nboth, nforw),
+                              nexog, ConstVector{ysteady});
         // form the shock realization
         ConstTwoDMatrix shocks_mat(nexog, nper, ConstVector{shocks});
         ConstTwoDMatrix vcov_mat(nexog, nexog, ConstVector{vcov});
diff --git a/dynare++/kord/faa_di_bruno.cc b/dynare++/kord/faa_di_bruno.cc
index 5cde7551b..3fa195b41 100644
--- a/dynare++/kord/faa_di_bruno.cc
+++ b/dynare++/kord/faa_di_bruno.cc
@@ -16,8 +16,8 @@ FaaDiBruno::calculate(const StackContainer<FGSTensor> &cont,
   out.zeros();
   for (int l = 1; l <= out.dimen(); l++)
     {
-      int mem_mb, p_size_mb;
-      int max = estimRefinment(out.getDims(), out.nrows(), l, mem_mb, p_size_mb);
+      int max, mem_mb, p_size_mb;
+      std::tie(max, mem_mb, p_size_mb) = estimRefinement(out.getDims(), out.nrows(), l);
       FoldedFineContainer fine_cont(cont, max);
       fine_cont.multAndAdd(l, f, out);
       JournalRecord recc(journal);
@@ -56,8 +56,8 @@ FaaDiBruno::calculate(const StackContainer<UGSTensor> &cont,
   out.zeros();
   for (int l = 1; l <= out.dimen(); l++)
     {
-      int mem_mb, p_size_mb;
-      int max = estimRefinment(out.getDims(), out.nrows(), l, mem_mb, p_size_mb);
+      int max, mem_mb, p_size_mb;
+      std::tie(max, mem_mb, p_size_mb) = estimRefinement(out.getDims(), out.nrows(), l);
       UnfoldedFineContainer fine_cont(cont, max);
       fine_cont.multAndAdd(l, f, out);
       JournalRecord recc(journal);
@@ -114,9 +114,8 @@ FaaDiBruno::calculate(const UnfoldedStackContainer &cont, const UGSContainer &g,
    If the right hand side is less than zero, we set |max| to 10, just to
    let it do something. */
 
-int
-FaaDiBruno::estimRefinment(const TensorDimens &tdims, int nr, int l,
-                           int &avmem_mb, int &tmpmem_mb)
+std::tuple<int, int, int>
+FaaDiBruno::estimRefinement(const TensorDimens &tdims, int nr, int l)
 {
   int nthreads = sthread::detach_thread_group::max_parallel_threads;
   long per_size1 = tdims.calcUnfoldMaxOffset();
@@ -142,7 +141,7 @@ FaaDiBruno::estimRefinment(const TensorDimens &tdims, int nr, int l,
         rec << " (decrease number of threads)";
       rec << endrec;
     }
-  avmem_mb = mem/1024/1024;
-  tmpmem_mb = nthreads*per_size/1024/1024;
-  return max;
+  int avmem_mb = mem/1024/1024;
+  int tmpmem_mb = nthreads*per_size/1024/1024;
+  return { max, avmem_mb, tmpmem_mb };
 }
diff --git a/dynare++/kord/faa_di_bruno.hh b/dynare++/kord/faa_di_bruno.hh
index aabcff027..15726eca2 100644
--- a/dynare++/kord/faa_di_bruno.hh
+++ b/dynare++/kord/faa_di_bruno.hh
@@ -17,6 +17,8 @@
 #include "sparse_tensor.hh"
 #include "gs_tensor.hh"
 
+#include <tuple>
+
 /* Nothing special here. See |@<|FaaDiBruno::calculate| folded sparse
    code@>| for reason of having |magic_mult|. */
 
@@ -37,7 +39,7 @@ public:
   void calculate(const UnfoldedStackContainer &cont, const UGSContainer &g,
                  UGSTensor &out);
 protected:
-  int estimRefinment(const TensorDimens &tdims, int nr, int l, int &avmem_mb, int &tmpmem_mb);
+  std::tuple<int,int,int> estimRefinement(const TensorDimens &tdims, int nr, int l);
   constexpr static double magic_mult = 1.5;
 };
 
diff --git a/dynare++/kord/kord_exception.hh b/dynare++/kord/kord_exception.hh
index 24857ee40..dd52cd853 100644
--- a/dynare++/kord/kord_exception.hh
+++ b/dynare++/kord/kord_exception.hh
@@ -4,12 +4,12 @@
 
 /* This is a simple code defining an exception and two convenience macros. */
 
+#include <string>
+#include <iostream>
+
 #ifndef KORD_EXCEPTION_H
 #define KORD_EXCEPTION_H
 
-#include <cstring>
-#include <cstdio>
-
 #define KORD_RAISE(mes)                         \
   throw KordException(__FILE__, __LINE__, mes);
 
@@ -25,31 +25,27 @@
 class KordException
 {
 protected:
-  char fname[50];
+  std::string fname;
   int lnum;
-  char message[500];
+  std::string message;
   int cd;
 public:
-  KordException(const char *f, int l, const char *mes, int c = 255)
+  KordException(std::string f, int l, std::string mes, int c = 255)
+    : fname{std::move(f)}, lnum{l}, message{std::move(mes)}, cd{c}
   {
-    strncpy(fname, f, 50); fname[49] = '\0';
-    strncpy(message, mes, 500); message[499] = '\0';
-    lnum = l;
-    cd = c;
   }
-  virtual ~KordException()
-  = default;
+  virtual ~KordException() = default;
   virtual void
   print() const
   {
-    printf("At %s:%d:(%d):%s\n", fname, lnum, cd, message);
+    std::cout << "At " << fname << ':' << lnum << ":(" << cd << "):" << message << '\n';
   }
   virtual int
   code() const
   {
     return cd;
   }
-  const char *
+  const std::string &
   get_message() const
   {
     return message;
@@ -57,8 +53,8 @@ public:
 };
 
 // |KordException| error code definitions
-#define KORD_FP_NOT_CONV 254
-#define KORD_FP_NOT_FINITE 253
-#define KORD_MD_NOT_STABLE 252
+constexpr int KORD_FP_NOT_CONV = 254;
+constexpr int KORD_FP_NOT_FINITE = 253;
+constexpr int KORD_MD_NOT_STABLE = 252;
 
 #endif
diff --git a/dynare++/kord/korder.hh b/dynare++/kord/korder.hh
index 953d9f5bd..9c7aee1ee 100644
--- a/dynare++/kord/korder.hh
+++ b/dynare++/kord/korder.hh
@@ -4,7 +4,7 @@
 
 /*
   The main purpose of this file is to implement a perturbation method
-  algorithm for an SDGE model for higher order approximations. The input
+  algorithm for an DSGE model for higher order approximations. The input
   of the algorithm are sparse tensors as derivatives of the dynamic
   system, then dimensions of vector variables, then the first order
   approximation to the decision rule and finally a covariance matrix of
@@ -76,10 +76,12 @@ public:
 
 /* The |PartitionY| class defines the partitioning of state variables
    $y$. The vector $y$, and subvector $y^*$, and $y^{**}$ are defined.
+
    $$y=\left[\matrix{\hbox{static}\cr\hbox{predeter}\cr\hbox{both}\cr
    \hbox{forward}}\right],\quad
    y^*=\left[\matrix{\hbox{predeter}\cr\hbox{both}}\right],\quad
    y^{**}=\left[\matrix{\hbox{both}\cr\hbox{forward}}\right],$$
+
    where ``static'' means variables appearing only at time $t$,
    ``predeter'' means variables appearing at time $t-1$, but not at
    $t+1$, ``both'' means variables appearing both at $t-1$ and $t+1$
@@ -180,60 +182,51 @@ public:
   }
 };
 
-/* Here we have the class for the higher order approximations. It
-   contains the following data:
+/* Here we have the class for the higher order approximations.
 
-   \halign{\kern\parindent\vrule height12pt width0pt
-   \vtop{\hsize=4cm\noindent\raggedright #}&\kern0.5cm\vtop{\hsize=10cm\noindent #}\cr
-   variable sizes ypart& |PartitionY| struct maintaining partitions of
-   $y$, see |@<|PartitionY| struct declaration@>|\cr
-   tensor variable dimension |nvs|& variable sizes of all tensors in
-   containers, sizes of $y^*$, $u$, $u'$@q'@> and $\sigma$\cr
-   tensor containers & folded and unfolded containers for $g$, $g_{y^*}$,
-   $g_{y^**}$ (the latter two collect appropriate subtensors of $g$, they
-   do not allocate any new space), $G$, $G$ stack, $Z$ stack\cr
-   dynamic model derivatives & just a reference to the container of
-   sparse tensors of the system derivatives, lives outside the class\cr
-   moments & both folded and unfolded normal moment containers, both are
-   calculated at initialization\cr
-   matrices & matrix $A$, matrix $S$, and matrix $B$, see |@<|MatrixA| class
-   declaration@>| and |@<|MatrixB| class declaration@>|\cr
-   }
+   It contains the following data:
 
-   \kern 0.4cm
+   - variable sizes ypart: |PartitionY| struct maintaining partitions of
+     $y$, see |@<|PartitionY| struct declaration@>|
+   - tensor variable dimension |nvs| : variable sizes of all tensors in
+     containers, sizes of $y^*$, $u$, $u'$ and $\sigma$
+   - tensor containers: folded and unfolded containers for $g$, $g_{y^*}$,
+     $g_{y^**}$ (the latter two collect appropriate subtensors of $g$, they
+     do not allocate any new space), $G$, $G$ stack, $Z$ stack
+   - dynamic model derivatives: just a reference to the container of
+     sparse tensors of the system derivatives, lives outside the class
+   - moments: both folded and unfolded normal moment containers, both are
+     calculated at initialization
+   - matrices: matrix $A$, matrix $S$, and matrix $B$, see |@<|MatrixA| class
+     declaration@>| and |@<|MatrixB| class declaration@>
 
    The methods are the following:
-   \halign{\kern\parindent\vrule height12pt width0pt
-   \vtop{\hsize=4cm\noindent\raggedright #}&\kern0.5cm\vtop{\hsize=10cm\noindent #}\cr
-   member access & we declare template methods for accessing containers
-   depending on |fold| and |unfold| flag, we implement their
-   specializations\cr
-   |performStep| & this performs $k$-order step provided that $k=2$ or
-   the $k-1$-th step has been run, this is the core method\cr
-   |check| & this calculates residuals of all solved equations for
-   $k$-order and reports their sizes, it is runnable after $k$-order
-   |performStep| has been run\cr
-   |insertDerivative| & inserts a $g$ derivative to the $g$ container and
-   also creates subtensors and insert them to $g_{y^*}$ and $g_{y^{**}}$
-   containers\cr
-   |sylvesterSolve| & solve the sylvester equation (templated fold, and
-   unfold)\cr
-   |faaDiBrunoZ| & calculates derivatives of $F$ by Faa Di Bruno for the
-   sparse container of system derivatives and $Z$ stack container\cr
-   |faaDiBrunoG| & calculates derivatives of $G$ by Faa Di Bruno for the
-   dense container $g^{**}$ and $G$ stack\cr
-   |recover_y| & recovers $g_{y^{*i}}$\cr
-   |recover_yu| & recovers $g_{y^{*i}u^j}$\cr
-   |recover_ys| & recovers $g_{y^{*i}\sigma^j}$\cr
-   |recover_yus| & recovers $g_{y^{*i}u^j\sigma^k}$\cr
-   |recover_s| & recovers $g_{\sigma^i}$\cr
-   |fillG| & calculates specified derivatives of $G$ and inserts them to
-   the container\cr
-   |calcE_ijk|& calculates $E_{ijk}$\cr
-   |calcD_ijk|& calculates $D_{ijk}$\cr
-   }
-
-   \kern 0.3cm
+   - member access: we declare template methods for accessing containers
+     depending on |fold| and |unfold| flag, we implement their
+     specializations
+   - |performStep|: this performs $k$-order step provided that $k=2$ or
+     the $k-1$-th step has been run, this is the core method
+   - |check|: this calculates residuals of all solved equations for
+     $k$-order and reports their sizes, it is runnable after $k$-order
+     |performStep| has been run
+   - |insertDerivative|: inserts a $g$ derivative to the $g$ container and
+     also creates subtensors and insert them to $g_{y^*}$ and $g_{y^{**}}$
+     containers
+   - |sylvesterSolve|: solve the sylvester equation (templated fold, and
+     unfold)
+   - |faaDiBrunoZ|: calculates derivatives of $F$ by Faa Di Bruno for the
+     sparse container of system derivatives and $Z$ stack container
+   - |faaDiBrunoG|: calculates derivatives of $G$ by Faa Di Bruno for the
+     dense container $g^{**}$ and $G$ stack
+   - |recover_y|: recovers $g_{y^{*i}}$
+   - |recover_yu|: recovers $g_{y^{*i}u^j}$
+   - |recover_ys|: recovers $g_{y^{*i}\sigma^j}$
+   - |recover_yus|: recovers $g_{y^{*i}u^j\sigma^k}$
+   - |recover_s|: recovers $g_{\sigma^i}$
+   - |fillG|: calculates specified derivatives of $G$ and inserts them to
+     the container
+   - |calcE_ijk|: calculates $E_{ijk}$
+   - |calcD_ijk|: calculates $D_{ijk}$
 
    Most of the code is templated, and template types are calculated in
    |ctraits|. So all templated methods get a template argument |T|, which
@@ -839,8 +832,7 @@ KOrder::check(int dim) const
       auto r = faaDiBrunoZ<t>(sym);
       double err = r->getData().getMax();
       JournalRecord(journal) << "\terror for symmetry " << sym << "\tis " << err << endrec;
-      if (err > maxerror)
-        maxerror = err;
+      maxerror = std::max(err, maxerror);
     }
 
   // check for $F_{y^iu^ju'^k}+D_{ijk}+E_{ijk}=0$
@@ -872,8 +864,7 @@ KOrder::check(int dim) const
   double err = r->getData().getMax();
   Symmetry sym{0, 0, 0, dim};
   JournalRecord(journal) << "\terror for symmetry " << sym << "\tis " << err << endrec;
-  if (err > maxerror)
-    maxerror = err;
+  maxerror = std::max(err, maxerror);
 
   return maxerror;
 }
diff --git a/dynare++/kord/normal_conjugate.cc b/dynare++/kord/normal_conjugate.cc
index a8e231b44..192ba88e0 100644
--- a/dynare++/kord/normal_conjugate.cc
+++ b/dynare++/kord/normal_conjugate.cc
@@ -29,10 +29,6 @@ NormalConj::NormalConj(const ConstTwoDMatrix &ydata)
     }
 }
 
-NormalConj::NormalConj(const NormalConj &nc)
-   
-= default;
-
 // |NormalConj::update| one observation code
 /* The method performs the following:
    $$\eqalign{
@@ -71,7 +67,7 @@ NormalConj::update(const ConstTwoDMatrix &ydata)
 void
 NormalConj::update(const NormalConj &nc)
 {
-  double wold = ((double) kappa)/(kappa+nc.kappa);
+  double wold = static_cast<double>(kappa)/(kappa+nc.kappa);
   double wnew = 1-wold;
 
   mu.mult(wold);
@@ -94,7 +90,7 @@ NormalConj::getVariance(TwoDMatrix &v) const
 {
   if (nu > getDim()+1)
     {
-      v = (const TwoDMatrix &) lambda;
+      v = const_cast<const TwoDMatrix &>(lambda);
       v.mult(1.0/(nu-getDim()-1));
     }
   else
diff --git a/dynare++/kord/normal_conjugate.hh b/dynare++/kord/normal_conjugate.hh
index d80d178c0..54637c36a 100644
--- a/dynare++/kord/normal_conjugate.hh
+++ b/dynare++/kord/normal_conjugate.hh
@@ -54,10 +54,10 @@ public:
      constructor. */
   NormalConj(int d);
   NormalConj(const ConstTwoDMatrix &ydata);
-  NormalConj(const NormalConj &nc);
+  NormalConj(const NormalConj &) = default;
+  NormalConj(NormalConj &&) = default;
 
-  virtual ~NormalConj()
-  = default;
+  virtual ~NormalConj() = default;
   void update(const ConstVector &y);
   void update(const ConstTwoDMatrix &ydata);
   void update(const NormalConj &nc);
diff --git a/dynare++/kord/tests.cc b/dynare++/kord/tests.cc
index ea1c9d5da..871aca6cf 100644
--- a/dynare++/kord/tests.cc
+++ b/dynare++/kord/tests.cc
@@ -309,8 +309,7 @@ TestRunnable::korder_unfold_fold(int maxdim, int unfold_dim,
                 << static_cast<double>(checktime)/CLOCKS_PER_SEC << '\n'
                 << "\tmax error in step dim=" << d << ":      " << std::setprecision(6) << err
                 << std::endl;
-      if (maxerror < err)
-        maxerror = err;
+      maxerror = std::max(err, maxerror);
     }
   // perform folded steps until maxdim
   if (unfold_dim < maxdim)
@@ -335,8 +334,7 @@ TestRunnable::korder_unfold_fold(int maxdim, int unfold_dim,
                     << static_cast<double>(checktime)/CLOCKS_PER_SEC << '\n'
                     << "\tmax error in step dim=" << d << ":      " << std::setprecision(6) << err
                     << std::endl;
-          if (maxerror < err)
-            maxerror = err;
+          maxerror = std::max(err, maxerror);
         }
     }
   return maxerror;
diff --git a/dynare++/src/dynare_exception.hh b/dynare++/src/dynare_exception.hh
index 6f0beb049..a0b985f2c 100644
--- a/dynare++/src/dynare_exception.hh
+++ b/dynare++/src/dynare_exception.hh
@@ -6,37 +6,26 @@
 #define DYNARE_EXCEPTION_H
 
 #include <string>
-#include <cstring>
+#include <utility>
 
 class DynareException
 {
-  char *mes;
+  std::string mes;
 public:
-  DynareException(const char *m, const char *fname, int line, int col)
+  DynareException(const std::string &m, const std::string &fname, int line, int col)
+    : mes{"Parse error at " + fname + ", line " + std::to_string(line) + ", column " +
+          std::to_string(col) + ": " + m}
   {
-    mes = new char[strlen(m) + strlen(fname) + 100];
-    sprintf(mes, "Parse error at %s, line %d, column %d: %s", fname, line, col, m);
   }
-  DynareException(const char *fname, int line, const std::string &m)
+  DynareException(const std::string &fname, int line, const std::string &m)
+    : mes{fname + ':' + std::to_string(line) + ": " + m}
   {
-    mes = new char[m.size() + strlen(fname) + 50];
-    sprintf(mes, "%s:%d: %s", fname, line, m.c_str());
   }
-  DynareException(const char *m, int offset)
+  DynareException(const std::string &m, int offset)
+    : mes{"Parse error in provided string at offset " + std::to_string(offset) + ": " + m}
   {
-    mes = new char[strlen(m) + 100];
-    sprintf(mes, "Parse error in provided string at offset %d: %s", offset, m);
   }
-  DynareException(const DynareException &e)
-    : mes(new char[strlen(e.mes)+1])
-  {
-    strcpy(mes, e.mes);
-  }
-  virtual ~DynareException()
-  {
-    delete [] mes;
-  }
-  const char *
+  const std::string &
   message() const
   {
     return mes;
diff --git a/dynare++/src/main.cc b/dynare++/src/main.cc
index f6b55552d..1dcb58c46 100644
--- a/dynare++/src/main.cc
+++ b/dynare++/src/main.cc
@@ -195,7 +195,7 @@ main(int argc, char **argv)
     }
   catch (const DynareException &e)
     {
-      printf("Caught Dynare exception: %s\n", e.message());
+      printf("Caught Dynare exception: %s\n", e.message().c_str());
       return 255;
     }
   catch (const ogu::Exception &e)
diff --git a/dynare++/tl/cc/stack_container.hh b/dynare++/tl/cc/stack_container.hh
index 86ccfad85..28e93401a 100644
--- a/dynare++/tl/cc/stack_container.hh
+++ b/dynare++/tl/cc/stack_container.hh
@@ -3,7 +3,7 @@
 // Stack of containers.
 
 /* Here we develop abstractions for stacked containers of tensors. For
-   instance, in perturbation methods for SDGE we need function
+   instance, in perturbation methods for DSGE we need function
    $$z(y,u,u',\sigma)=\left[\matrix{G(y,u,u',\sigma)\cr g(y,u,\sigma)\cr y\cr u}\right]$$
    and we need to calculate one step of Faa Di Bruno formula
    $$\left[B_{s^k}\right]_{\alpha_1\ldots\alpha_l}=\left[f_{z^l}\right]_{\beta_1\ldots\beta_l}
@@ -319,7 +319,7 @@ protected:
 };
 
 /* Here is the specialization of the |StackContainer|. We implement
-   here the $z$ needed in SDGE context. We implement |getType| and define
+   here the $z$ needed in DSGE context. We implement |getType| and define
    a constructor feeding the data and sizes.
 
    Note that it has two containers, the first is dependent on four
@@ -401,7 +401,7 @@ public:
 };
 
 /* Here we have another specialization of container used in context of
-   SDGE. We define a container for
+   DSGE. We define a container for
    $$G(y,u,u',\sigma)=g^{**}(g^*(y,u,\sigma),u',\sigma)$$
 
    For some reason, the symmetry of $g^{**}$ has length $4$ although it
diff --git a/mex/sources/k_order_perturbation/dynamic_abstract_class.cc b/mex/sources/k_order_perturbation/dynamic_abstract_class.cc
index d4fd0f72f..76bcb6461 100644
--- a/mex/sources/k_order_perturbation/dynamic_abstract_class.cc
+++ b/mex/sources/k_order_perturbation/dynamic_abstract_class.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2010-2014 Dynare Team
+ * Copyright (C) 2010-2019 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -50,7 +50,7 @@ DynamicModelAC::unpackSparseMatrixAndCopyIntoTwoDMatData(mxArray *sparseMat, Two
   int output_row = 0;
 
   for (int i = 0; i < totalCols; i++)
-    for (int j = 0; j < (int) (colIdxVector[i+1]-colIdxVector[i]); j++, rind++)
+    for (int j = 0; j < static_cast<int>((colIdxVector[i+1]-colIdxVector[i])); j++, rind++)
       {
         tdm->get(output_row, 0) = rowIdxVector[rind] + 1;
         tdm->get(output_row, 1) = i + 1;
@@ -62,7 +62,7 @@ DynamicModelAC::unpackSparseMatrixAndCopyIntoTwoDMatData(mxArray *sparseMat, Two
      derivative is symbolically not zero but numerically zero at the evaluation
      point), then fill in the matrix with empty entries, that will be
      recognized as such by KordpDynare::populateDerivativesContainer() */
-  while (output_row < (int) sizeRowIdxVector)
+  while (output_row < static_cast<int>(sizeRowIdxVector))
     {
       tdm->get(output_row, 0) = 0;
       tdm->get(output_row, 1) = 0;
diff --git a/mex/sources/k_order_perturbation/dynamic_abstract_class.hh b/mex/sources/k_order_perturbation/dynamic_abstract_class.hh
index e7f371d82..9cdd2621c 100644
--- a/mex/sources/k_order_perturbation/dynamic_abstract_class.hh
+++ b/mex/sources/k_order_perturbation/dynamic_abstract_class.hh
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2010-2014 Dynare Team
+ * Copyright (C) 2010-2019 Dynare Team
  *
  * This file is part of Dynare.
  *
diff --git a/mex/sources/k_order_perturbation/dynamic_dll.cc b/mex/sources/k_order_perturbation/dynamic_dll.cc
index 60904e82c..719a93a46 100644
--- a/mex/sources/k_order_perturbation/dynamic_dll.cc
+++ b/mex/sources/k_order_perturbation/dynamic_dll.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2008-2012 Dynare Team
+ * Copyright (C) 2008-2019 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -19,8 +19,6 @@
 
 #include "dynamic_dll.hh"
 
-#include <sstream>
-
 DynamicModelDLL::DynamicModelDLL(const std::string &modName) noexcept(false)
 {
   std::string fName;
@@ -33,22 +31,21 @@ DynamicModelDLL::DynamicModelDLL(const std::string &modName) noexcept(false)
     {
 #if defined(__CYGWIN32__) || defined(_WIN32)
       dynamicHinstance = LoadLibrary(fName.c_str());
-      if (dynamicHinstance == nullptr)
+      if (!dynamicHinstance)
         throw 1;
-      ntt = (int *) GetProcAddress(dynamicHinstance, "ntt");
-      dynamic_resid_tt = (dynamic_tt_fct) GetProcAddress(dynamicHinstance, "dynamic_resid_tt");
-      dynamic_resid = (dynamic_resid_fct) GetProcAddress(dynamicHinstance, "dynamic_resid");
-      dynamic_g1_tt = (dynamic_tt_fct) GetProcAddress(dynamicHinstance, "dynamic_g1_tt");
-      dynamic_g1 = (dynamic_g1_fct) GetProcAddress(dynamicHinstance, "dynamic_g1");
-      dynamic_g2_tt = (dynamic_tt_fct) GetProcAddress(dynamicHinstance, "dynamic_g2_tt");
-      dynamic_g2 = (dynamic_g2_fct) GetProcAddress(dynamicHinstance, "dynamic_g2");
-      dynamic_g3_tt = (dynamic_tt_fct) GetProcAddress(dynamicHinstance, "dynamic_g3_tt");
-      dynamic_g3 = (dynamic_g3_fct) GetProcAddress(dynamicHinstance, "dynamic_g3");
-      if (ntt == nullptr
-          || dynamic_resid_tt == nullptr || dynamic_resid == nullptr
-          || dynamic_g1_tt == nullptr || dynamic_g1 == nullptr
-          || dynamic_g2_tt == nullptr || dynamic_g2 == nullptr
-          || dynamic_g3_tt == nullptr || dynamic_g3 == nullptr)
+      ntt = reinterpret_cast<int *>(GetProcAddress(dynamicHinstance, "ntt"));
+      dynamic_resid_tt = reinterpret_cast<dynamic_tt_fct>(GetProcAddress(dynamicHinstance, "dynamic_resid_tt"));
+      dynamic_resid = reinterpret_cast<dynamic_resid_fct>(GetProcAddress(dynamicHinstance, "dynamic_resid"));
+      dynamic_g1_tt = reinterpret_cast<dynamic_tt_fct>(GetProcAddress(dynamicHinstance, "dynamic_g1_tt"));
+      dynamic_g1 = reinterpret_cast<dynamic_g1_fct>(GetProcAddress(dynamicHinstance, "dynamic_g1"));
+      dynamic_g2_tt = reinterpret_cast<dynamic_tt_fct>(GetProcAddress(dynamicHinstance, "dynamic_g2_tt"));
+      dynamic_g2 = reinterpret_cast<dynamic_g2_fct>(GetProcAddress(dynamicHinstance, "dynamic_g2"));
+      dynamic_g3_tt = reinterpret_cast<dynamic_tt_fct>(GetProcAddress(dynamicHinstance, "dynamic_g3_tt"));
+      dynamic_g3 = reinterpret_cast<dynamic_g3_fct>(GetProcAddress(dynamicHinstance, "dynamic_g3"));
+      if (!ntt || !dynamic_resid_tt || !dynamic_resid
+          || !dynamic_g1_tt || !dynamic_g1
+          || !dynamic_g2_tt || !dynamic_g2
+          || !dynamic_g3_tt || !dynamic_g3)
         {
           FreeLibrary(dynamicHinstance); // Free the library
           throw 2;
@@ -60,20 +57,19 @@ DynamicModelDLL::DynamicModelDLL(const std::string &modName) noexcept(false)
           std::cerr << dlerror() << std::endl;
           throw 1;
         }
-      ntt = (int *) dlsym(dynamicHinstance, "ntt");
-      dynamic_resid_tt = (dynamic_tt_fct) dlsym(dynamicHinstance, "dynamic_resid_tt");
-      dynamic_resid = (dynamic_resid_fct) dlsym(dynamicHinstance, "dynamic_resid");
-      dynamic_g1_tt = (dynamic_tt_fct) dlsym(dynamicHinstance, "dynamic_g1_tt");
-      dynamic_g1 = (dynamic_g1_fct) dlsym(dynamicHinstance, "dynamic_g1");
-      dynamic_g2_tt = (dynamic_tt_fct) dlsym(dynamicHinstance, "dynamic_g2_tt");
-      dynamic_g2 = (dynamic_g2_fct) dlsym(dynamicHinstance, "dynamic_g2");
-      dynamic_g3_tt = (dynamic_tt_fct) dlsym(dynamicHinstance, "dynamic_g3_tt");
-      dynamic_g3 = (dynamic_g3_fct) dlsym(dynamicHinstance, "dynamic_g3");
-      if (ntt == nullptr
-          || dynamic_resid_tt == nullptr || dynamic_resid == nullptr
-          || dynamic_g1_tt == nullptr || dynamic_g1 == nullptr
-          || dynamic_g2_tt == nullptr || dynamic_g2 == nullptr
-          || dynamic_g3_tt == nullptr || dynamic_g3 == nullptr)
+      ntt = reinterpret_cast<int *>(dlsym(dynamicHinstance, "ntt"));
+      dynamic_resid_tt = reinterpret_cast<dynamic_tt_fct>(dlsym(dynamicHinstance, "dynamic_resid_tt"));
+      dynamic_resid = reinterpret_cast<dynamic_resid_fct>(dlsym(dynamicHinstance, "dynamic_resid"));
+      dynamic_g1_tt = reinterpret_cast<dynamic_tt_fct>(dlsym(dynamicHinstance, "dynamic_g1_tt"));
+      dynamic_g1 = reinterpret_cast<dynamic_g1_fct>(dlsym(dynamicHinstance, "dynamic_g1"));
+      dynamic_g2_tt = reinterpret_cast<dynamic_tt_fct>(dlsym(dynamicHinstance, "dynamic_g2_tt"));
+      dynamic_g2 = reinterpret_cast<dynamic_g2_fct>(dlsym(dynamicHinstance, "dynamic_g2"));
+      dynamic_g3_tt = reinterpret_cast<dynamic_tt_fct>(dlsym(dynamicHinstance, "dynamic_g3_tt"));
+      dynamic_g3 = reinterpret_cast<dynamic_g3_fct>(dlsym(dynamicHinstance, "dynamic_g3"));
+      if (!ntt || !dynamic_resid_tt || !dynamic_resid
+          || !dynamic_g1_tt || !dynamic_g1
+          || !dynamic_g2_tt || !dynamic_g2
+          || !dynamic_g3_tt || !dynamic_g3)
         {
           dlclose(dynamicHinstance); // Free the library
           std::cerr << dlerror() << std::endl;
@@ -84,14 +80,13 @@ DynamicModelDLL::DynamicModelDLL(const std::string &modName) noexcept(false)
     }
   catch (int i)
     {
-      std::ostringstream msg;
-      msg << "Error when loading " << fName << " (";
+      std::string msg{"Error when loading " + fName + " ("};
       if (i == 1)
-        msg << "can't dynamically load the file";
+        msg += "can't dynamically load the file";
       if (i == 2)
-        msg << "can't locate the relevant dynamic symbols within the MEX file";
-      msg << ")";
-      throw DynareException(__FILE__, __LINE__, msg.str());
+        msg += "can't locate the relevant dynamic symbols within the MEX file";
+      msg += ")";
+      throw DynareException(__FILE__, __LINE__, msg);
     }
   catch (...)
     {
@@ -117,7 +112,7 @@ void
 DynamicModelDLL::eval(const Vector &y, const Vector &x, const Vector &modParams, const Vector &ySteady,
                       Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3) noexcept(false)
 {
-  double *T = (double *) malloc(sizeof(double) * (*ntt));
+  double *T = new double[*ntt];
   dynamic_resid_tt(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, T);
   dynamic_resid(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, T, residual.base());
   if (g1 || g2 || g3)
@@ -133,5 +128,5 @@ DynamicModelDLL::eval(const Vector &y, const Vector &x, const Vector &modParams,
       dynamic_g3_tt(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, T);
       dynamic_g3(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, T, g3->base());
     }
-  free(T);
+  delete[] T;
 }
diff --git a/mex/sources/k_order_perturbation/dynamic_dll.hh b/mex/sources/k_order_perturbation/dynamic_dll.hh
index 9385afa35..3e0568029 100644
--- a/mex/sources/k_order_perturbation/dynamic_dll.hh
+++ b/mex/sources/k_order_perturbation/dynamic_dll.hh
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2008-2017 Dynare Team
+ * Copyright (C) 2008-2019 Dynare Team
  *
  * This file is part of Dynare.
  *
diff --git a/mex/sources/k_order_perturbation/dynamic_m.cc b/mex/sources/k_order_perturbation/dynamic_m.cc
index 6e37c3f54..1ef2dfa8d 100644
--- a/mex/sources/k_order_perturbation/dynamic_m.cc
+++ b/mex/sources/k_order_perturbation/dynamic_m.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2010-2011 Dynare Team
+ * Copyright (C) 2010-2019 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -17,6 +17,8 @@
  * along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
  */
 
+#include <algorithm>
+
 #include "dynamic_m.hh"
 
 DynamicModelMFile::DynamicModelMFile(const std::string &modName) noexcept(false) :
@@ -36,20 +38,21 @@ DynamicModelMFile::eval(const Vector &y, const Vector &x, const Vector &modParam
   prhs[3] = mxCreateDoubleMatrix(ySteady.length(), 1, mxREAL);
   prhs[4] = mxCreateDoubleScalar(1.0);
 
-  memcpy(mxGetData(prhs[0]), (void *) y.base(), y.length()*sizeof(double));
-  memcpy(mxGetData(prhs[1]), (void *) x.base(), x.length()*sizeof(double));
-  memcpy(mxGetData(prhs[2]), (void *) modParams.base(), modParams.length()*sizeof(double));
-  memcpy(mxGetData(prhs[3]), (void *) ySteady.base(), ySteady.length()*sizeof(double));
+  std::copy_n(y.base(), y.length(), mxGetPr(prhs[0]));
+  std::copy_n(x.base(), x.length(), mxGetPr(prhs[1]));
+  std::copy_n(modParams.base(), modParams.length(), mxGetPr(prhs[2]));
+  std::copy_n(ySteady.base(), ySteady.length(), mxGetPr(prhs[3]));
 
   int retVal = mexCallMATLAB(nlhs_dynamic, plhs, nrhs_dynamic, prhs, DynamicMFilename.c_str());
   if (retVal != 0)
     throw DynareException(__FILE__, __LINE__, "Trouble calling " + DynamicMFilename);
 
   residual = Vector{plhs[0]};
-  copyDoubleIntoTwoDMatData(mxGetPr(plhs[1]), g1, (int) mxGetM(plhs[1]), (int) mxGetN(plhs[1]));
-  if (g2 != nullptr)
+  copyDoubleIntoTwoDMatData(mxGetPr(plhs[1]), g1, static_cast<int>(mxGetM(plhs[1])),
+                            static_cast<int>(mxGetN(plhs[1])));
+  if (g2)
     unpackSparseMatrixAndCopyIntoTwoDMatData(plhs[2], g2);
-  if (g3 != nullptr)
+  if (g3)
     unpackSparseMatrixAndCopyIntoTwoDMatData(plhs[3], g3);
 
   for (int i = 0; i < nrhs_dynamic; i++)
diff --git a/mex/sources/k_order_perturbation/dynamic_m.hh b/mex/sources/k_order_perturbation/dynamic_m.hh
index 48c2b396b..4632f3a33 100644
--- a/mex/sources/k_order_perturbation/dynamic_m.hh
+++ b/mex/sources/k_order_perturbation/dynamic_m.hh
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2010-2017 Dynare Team
+ * Copyright (C) 2010-2019 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -33,8 +33,8 @@ class DynamicModelMFile : public DynamicModelAC
 {
 private:
   const std::string DynamicMFilename;
-  const static int nlhs_dynamic = 4;
-  const static int nrhs_dynamic = 5;
+  constexpr static int nlhs_dynamic = 4;
+  constexpr static int nrhs_dynamic = 5;
 public:
   explicit DynamicModelMFile(const std::string &modName) noexcept(false);
   virtual ~DynamicModelMFile() = default;
diff --git a/mex/sources/k_order_perturbation/k_ord_dynare.cc b/mex/sources/k_order_perturbation/k_ord_dynare.cc
index 00e0c2d78..be1dd34ff 100644
--- a/mex/sources/k_order_perturbation/k_ord_dynare.cc
+++ b/mex/sources/k_order_perturbation/k_ord_dynare.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2008-2017 Dynare Team
+ * Copyright (C) 2008-2019 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -44,8 +44,8 @@ KordpDynare::KordpDynare(const std::vector<std::string> &endo, int num_endo,
   nYs{npred + nboth}, nYss{nboth + nforw}, nY{num_endo}, nJcols{jcols}, NNZD{nnzd}, nSteps{nsteps},
   nOrder{norder}, journal{jr}, ySteady{ysteady}, params{inParams}, vCov{vcov},
   md{1}, dnl{*this, endo}, denl{*this, exo}, dsnl{*this, dnl, denl}, ss_tol{sstol}, varOrder{var_order},
-  ll_Incidence{llincidence}, qz_criterium{criterium}, g1p{nullptr},
-  g2p{nullptr}, g3p{nullptr}, dynamicModelFile{move(dynamicModelFile_arg)}
+  ll_Incidence{llincidence}, qz_criterium{criterium},
+  dynamicModelFile{move(dynamicModelFile_arg)}
 {
   ReorderDynareJacobianIndices();
 
@@ -61,13 +61,15 @@ KordpDynare::KordpDynare(const std::vector<std::string> &endo, int num_endo,
                          const int nsteps, int norder,
                          Journal &jr, std::unique_ptr<DynamicModelAC> dynamicModelFile_arg, double sstol,
                          const std::vector<int> &var_order, const TwoDMatrix &llincidence, double criterium,
-                         TwoDMatrix *g1_arg, TwoDMatrix *g2_arg, TwoDMatrix *g3_arg) noexcept(false) :
+                         std::unique_ptr<TwoDMatrix> g1_arg, std::unique_ptr<TwoDMatrix> g2_arg,
+                         std::unique_ptr<TwoDMatrix> g3_arg) noexcept(false) :
   nStat{nstat}, nBoth{nboth}, nPred{npred}, nForw{nforw}, nExog{nexog}, nPar{npar},
   nYs{npred + nboth}, nYss{nboth + nforw}, nY{num_endo}, nJcols{jcols}, NNZD{nnzd}, nSteps{nsteps},
   nOrder{norder}, journal{jr}, ySteady{ysteady}, params{inParams}, vCov{vcov},
   md{1}, dnl{*this, endo}, denl{*this, exo}, dsnl{*this, dnl, denl}, ss_tol{sstol}, varOrder{var_order},
   ll_Incidence{llincidence}, qz_criterium{criterium},
-  g1p{g1_arg}, g2p{g2_arg}, g3p{g3_arg}, dynamicModelFile{move(dynamicModelFile_arg)}
+  g1p{std::move(g1_arg)}, g2p{std::move(g2_arg)}, g3p{std::move(g3_arg)},
+  dynamicModelFile{move(dynamicModelFile_arg)}
 {
   ReorderDynareJacobianIndices();
 
@@ -76,11 +78,6 @@ KordpDynare::KordpDynare(const std::vector<std::string> &endo, int num_endo,
     md.insert(std::make_unique<FSSparseTensor>(iord, nY+nYs+nYss+nExog, nY));
 }
 
-KordpDynare::~KordpDynare()
-{
-  // No need to manually delete tensors in "md", they are deleted by the TensorContainer destructor
-}
-
 void
 KordpDynare::solveDeterministicSteady()
 {
@@ -110,21 +107,21 @@ KordpDynare::evaluateSystem(Vector &out, const ConstVector &yym, const ConstVect
 void
 KordpDynare::calcDerivativesAtSteady()
 {
-  if (g1p == nullptr)
+  if (!g1p)
     {
-      g1p = new TwoDMatrix(nY, nJcols);
+      g1p = std::make_unique<TwoDMatrix>(nY, nJcols);
       g1p->zeros();
 
       if (nOrder > 1)
         {
           // allocate space for sparse Hessian
-          g2p = new TwoDMatrix((int) NNZD[1], 3);
+          g2p = std::make_unique<TwoDMatrix>(static_cast<int>(NNZD[1]), 3);
           g2p->zeros();
         }
 
       if (nOrder > 2)
         {
-          g3p = new TwoDMatrix((int) NNZD[2], 3);
+          g3p = std::make_unique<TwoDMatrix>(static_cast<int>(NNZD[2]), 3);
           g3p->zeros();
         }
 
@@ -136,23 +133,16 @@ KordpDynare::calcDerivativesAtSteady()
       Vector llxSteady(nJcols-nExog);
       LLxSteady(ySteady, llxSteady);
 
-      dynamicModelFile->eval(llxSteady, xx, params, ySteady, out, g1p, g2p, g3p);
-
+      dynamicModelFile->eval(llxSteady, xx, params, ySteady, out, g1p.get(), g2p.get(), g3p.get());
     }
 
   populateDerivativesContainer(*g1p, 1, JacobianIndices);
-  delete g1p;
 
   if (nOrder > 1)
-    {
-      populateDerivativesContainer(*g2p, 2, JacobianIndices);
-      delete g2p;
-    }
+    populateDerivativesContainer(*g2p, 2, JacobianIndices);
+
   if (nOrder > 2)
-    {
-      populateDerivativesContainer(*g3p, 3, JacobianIndices);
-      delete g3p;
-    }
+    populateDerivativesContainer(*g3p, 3, JacobianIndices);
 }
 
 /*******************************************************************************
@@ -167,22 +157,20 @@ KordpDynare::populateDerivativesContainer(const TwoDMatrix &g, int ord, const st
   IntSequence s(ord, 0);
 
   if (ord == 1)
-    {
-      for (int i = 0; i < g.ncols(); i++)
-        {
-          for (int j = 0; j < g.nrows(); j++)
-            {
-              double x;
-              if (s[0] < nJcols-nExog)
-                x = g.get(j, vOrder[s[0]]);
-              else
-                x = g.get(j, s[0]);
-              if (x != 0.0)
-                mdTi->insert(s, j, x);
-            }
-          s[0]++;
-        }
-    }
+    for (int i = 0; i < g.ncols(); i++)
+      {
+        for (int j = 0; j < g.nrows(); j++)
+          {
+            double x;
+            if (s[0] < nJcols-nExog)
+              x = g.get(j, vOrder[s[0]]);
+            else
+              x = g.get(j, s[0]);
+            if (x != 0.0)
+              mdTi->insert(s, j, x);
+          }
+        s[0]++;
+      }
   else if (ord == 2)
     {
       int nJcols1 = nJcols-nExog;
@@ -191,8 +179,8 @@ KordpDynare::populateDerivativesContainer(const TwoDMatrix &g, int ord, const st
         revOrder[vOrder[i]] = i;
       for (int i = 0; i < g.nrows(); i++)
         {
-          int j = (int) g.get(i, 0)-1; // hessian indices start with 1
-          int i1 = (int) g.get(i, 1) -1;
+          int j = static_cast<int>(g.get(i, 0))-1; // hessian indices start with 1
+          int i1 = static_cast<int>(g.get(i, 1)) -1;
           if (j < 0 || i1 < 0)
             continue; // Discard empty entries (see comment in DynamicModelAC::unpackSparseMatrix())
           int s0 = i1 / nJcols;
@@ -221,8 +209,8 @@ KordpDynare::populateDerivativesContainer(const TwoDMatrix &g, int ord, const st
         revOrder[vOrder[i]] = i;
       for (int i = 0; i < g.nrows(); i++)
         {
-          int j = (int) g.get(i, 0)-1;
-          int i1 = (int) g.get(i, 1) -1;
+          int j = static_cast<int>(g.get(i, 0))-1;
+          int i1 = static_cast<int>(g.get(i, 1)) -1;
           if (j < 0 || i1 < 0)
             continue; // Discard empty entries (see comment in DynamicModelAC::unpackSparseMatrix())
           int s0 = i1 / nJcols2;
@@ -269,14 +257,12 @@ KordpDynare::LLxSteady(const Vector &yS, Vector &llxSteady) noexcept(false)
     throw DynareException(__FILE__, __LINE__, "llxSteady has wrong size");
 
   for (int ll_row = 0; ll_row < ll_Incidence.nrows(); ll_row++)
-    {
-      // populate (non-sparse) vector with ysteady values
-      for (int i = 0; i < nY; i++)
-        {
-          if (ll_Incidence.get(ll_row, i))
-            llxSteady[((int) ll_Incidence.get(ll_row, i))-1] = yS[i];
-        }
-    }
+    // populate (non-sparse) vector with ysteady values
+    for (int i = 0; i < nY; i++)
+      {
+        if (ll_Incidence.get(ll_row, i))
+          llxSteady[static_cast<int>(ll_Incidence.get(ll_row, i))-1] = yS[i];
+      }
 }
 
 /************************************
@@ -318,7 +304,7 @@ KordpDynare::ReorderDynareJacobianIndices() noexcept(false)
       // reorder in orde-var order & populate temporary nEndo (sparse) vector with
       // the lag, current and lead blocks of the jacobian respectively
       for (i = 0; i < nY; i++)
-        tmp[i] = ((int) ll_Incidence.get(ll_row, varOrder[i]-1));
+        tmp[i] = static_cast<int>(ll_Incidence.get(ll_row, varOrder[i]-1));
       // write the reordered blocks back to the jacobian
       // in reverse order
       for (j = nY-1; j >= 0; j--)
@@ -340,7 +326,8 @@ KordpDynare::ReorderDynareJacobianIndices() noexcept(false)
 /*       DynareNameList class                                                         */
 /**************************************************************************************/
 
-DynareNameList::DynareNameList(const KordpDynare &dynare, std::vector<std::string> names_arg) : names(std::move(names_arg))
+DynareNameList::DynareNameList(const KordpDynare &dynare, std::vector<std::string> names_arg)
+  : names(std::move(names_arg))
 {
 }
 
diff --git a/mex/sources/k_order_perturbation/k_ord_dynare.hh b/mex/sources/k_order_perturbation/k_ord_dynare.hh
index f8325083c..4654f6129 100644
--- a/mex/sources/k_order_perturbation/k_ord_dynare.hh
+++ b/mex/sources/k_order_perturbation/k_ord_dynare.hh
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2008-2017 Dynare Team
+ * Copyright (C) 2008-2019 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -19,8 +19,10 @@
 
 #ifndef K_ORD_DYNARE3_H
 #define K_ORD_DYNARE3_H
+
 #include <vector>
 #include <memory>
+
 #include "t_container.hh"
 #include "sparse_tensor.hh"
 #include "decision_rule.hh"
@@ -115,9 +117,7 @@ class KordpDynare : public DynamicModel
   double qz_criterium;
   std::vector<int> JacobianIndices;
 
-  TwoDMatrix *g1p;
-  TwoDMatrix *g2p;
-  TwoDMatrix *g3p;
+  std::unique_ptr<TwoDMatrix> g1p, g2p, g3p;
 public:
   KordpDynare(const std::vector<std::string> &endo, int num_endo,
               const std::vector<std::string> &exo, int num_exo, int num_par,
@@ -134,9 +134,9 @@ public:
               const int nSteps, const int ord,
               Journal &jr, std::unique_ptr<DynamicModelAC> dynamicModelFile_arg, double sstol,
               const std::vector<int> &varOrder, const TwoDMatrix &ll_Incidence,
-              double qz_criterium, TwoDMatrix *g1_arg, TwoDMatrix *g2_arg, TwoDMatrix *g3_arg) noexcept(false);
+              double qz_criterium, std::unique_ptr<TwoDMatrix> g1_arg,
+              std::unique_ptr<TwoDMatrix> g2_arg, std::unique_ptr<TwoDMatrix> g3_arg) noexcept(false);
 
-  virtual ~KordpDynare();
   int
   nstat() const
   {
diff --git a/mex/sources/k_order_perturbation/k_order_perturbation.cc b/mex/sources/k_order_perturbation/k_order_perturbation.cc
index 696025194..233cf0798 100644
--- a/mex/sources/k_order_perturbation/k_order_perturbation.cc
+++ b/mex/sources/k_order_perturbation/k_order_perturbation.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2008-2017 Dynare Team
+ * Copyright (C) 2008-2019 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -34,14 +34,11 @@
 #include "dynamic_m.hh"
 #include "dynamic_dll.hh"
 
+#include <algorithm>
 #include <cmath>
-#include <cstring>
-#include <cctype>
 #include <cassert>
 
-#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)  // exclude mexFunction for other applications
-
-# include "dynmex.h"
+#include "dynmex.h"
 
 //////////////////////////////////////////////////////
 // Convert MATLAB Dynare endo and exo names array to a vector<string> array of string pointers
@@ -69,7 +66,7 @@ copy_derivatives(mxArray *destin, const Symmetry &sym, const FGSContainer &deriv
   int n = x.numRows();
   int m = x.numCols();
   mxArray *tmp = mxCreateDoubleMatrix(n, m, mxREAL);
-  memcpy(mxGetPr(tmp), x.getData().base(), n*m*sizeof(double));
+  std::copy_n(x.getData().base(), n*m, mxGetPr(tmp));
   mxSetField(destin, 0, fieldname.c_str(), tmp);
 }
 
@@ -85,7 +82,7 @@ extern "C" {
     const mxArray *dr = prhs[0];
     const mxArray *M_ = prhs[1];
     const mxArray *options_ = prhs[2];
-    int use_dll = (int) mxGetScalar(mxGetField(options_, 0, "use_dll"));
+    bool use_dll = mxGetScalar(mxGetField(options_, 0, "use_dll")) != 0;
 
     mxArray *mFname = mxGetField(M_, 0, "fname");
     if (!mxIsChar(mFname))
@@ -96,19 +93,14 @@ extern "C" {
     int kOrder;
     mxArray *mxFldp = mxGetField(options_, 0, "order");
     if (mxIsNumeric(mxFldp))
-      kOrder = (int) mxGetScalar(mxFldp);
+      kOrder = static_cast<int>(mxGetScalar(mxFldp));
     else
       kOrder = 1;
 
-    //if (kOrder == 1 && nlhs != 2)
-    //  DYN_MEX_FUNC_ERR_MSG_TXT("k_order_perturbation at order 1 requires exactly 2 arguments in output");
-    //else if (kOrder > 1 && nlhs != kOrder+2)
-    //  DYN_MEX_FUNC_ERR_MSG_TXT("k_order_perturbation at order > 1 requires exactly order+2 arguments in output");
-
     double qz_criterium = 1+1e-6;
     mxFldp = mxGetField(options_, 0, "qz_criterium");
     if (mxGetNumberOfElements(mxFldp) > 0 && mxIsNumeric(mxFldp))
-      qz_criterium = (double) mxGetScalar(mxFldp);
+      qz_criterium = mxGetScalar(mxFldp);
 
     mxFldp = mxGetField(M_, 0, "params");
     Vector modParams{mxFldp};
@@ -127,47 +119,45 @@ extern "C" {
       DYN_MEX_FUNC_ERR_MSG_TXT("The steady state vector contains NaN or Inf");
 
     mxFldp = mxGetField(M_, 0, "nstatic");
-    const int nStat = (int) mxGetScalar(mxFldp);
+    const int nStat = static_cast<int>(mxGetScalar(mxFldp));
     mxFldp = mxGetField(M_, 0, "npred");
-    const int nPred = (int) mxGetScalar(mxFldp);
+    const int nPred = static_cast<int>(mxGetScalar(mxFldp));
     mxFldp = mxGetField(M_, 0, "nspred");
-    const int nsPred = (int) mxGetScalar(mxFldp);
+    const int nsPred = static_cast<int>(mxGetScalar(mxFldp));
     mxFldp = mxGetField(M_, 0, "nboth");
-    const int nBoth = (int) mxGetScalar(mxFldp);
+    const int nBoth = static_cast<int>(mxGetScalar(mxFldp));
     mxFldp = mxGetField(M_, 0, "nfwrd");
-    const int nForw = (int) mxGetScalar(mxFldp);
+    const int nForw = static_cast<int>(mxGetScalar(mxFldp));
     mxFldp = mxGetField(M_, 0, "nsfwrd");
-    const int nsForw = (int) mxGetScalar(mxFldp);
+    const int nsForw = static_cast<int>(mxGetScalar(mxFldp));
 
     mxFldp = mxGetField(M_, 0, "exo_nbr");
-    const int nExog = (int) mxGetScalar(mxFldp);
+    const int nExog = static_cast<int>(mxGetScalar(mxFldp));
     mxFldp = mxGetField(M_, 0, "endo_nbr");
-    const int nEndo = (int) mxGetScalar(mxFldp);
+    const int nEndo = static_cast<int>(mxGetScalar(mxFldp));
     mxFldp = mxGetField(M_, 0, "param_nbr");
-    const int nPar = (int) mxGetScalar(mxFldp);
+    const int nPar = static_cast<int>(mxGetScalar(mxFldp));
 
     mxFldp = mxGetField(dr, 0, "order_var");
     auto dparams = mxGetPr(mxFldp);
-    npar = (int) mxGetM(mxFldp);
+    npar = static_cast<int>(mxGetM(mxFldp));
     if (npar != nEndo)
       DYN_MEX_FUNC_ERR_MSG_TXT("Incorrect number of input var_order vars.");
 
     std::vector<int> var_order_vp(nEndo);
     for (int v = 0; v < nEndo; v++)
-      var_order_vp[v] = (int) (*(dparams++));
+      var_order_vp[v] = static_cast<int>(*(dparams++));
 
     // the lag, current and lead blocks of the jacobian respectively
     mxFldp = mxGetField(M_, 0, "lead_lag_incidence");
-    npar = (int) mxGetN(mxFldp);
-    int nrows = (int) mxGetM(mxFldp);
+    npar = static_cast<int>(mxGetN(mxFldp));
+    int nrows = static_cast<int>(mxGetM(mxFldp));
 
     TwoDMatrix llincidence(nrows, npar, Vector{mxFldp});
     if (npar != nEndo)
-      {
-        std::ostringstream strstrm;
-        strstrm << "dynare:k_order_perturbation " << "Incorrect length of lead lag incidences: ncol=" << npar << " != nEndo=" << nEndo;
-        DYN_MEX_FUNC_ERR_MSG_TXT(strstrm.str().c_str());
-      }
+      DYN_MEX_FUNC_ERR_MSG_TXT("dynare:k_order_perturbation: Incorrect length of lead lag incidences: ncol="
+                               + std::to_string(npar) + " != nEndo=" + std::to_string(nEndo));
+
     //get NNZH =NNZD(2) = the total number of non-zero Hessian elements
     mxFldp = mxGetField(M_, 0, "NNZDerivatives");
     Vector NNZD{mxFldp};
@@ -177,42 +167,40 @@ extern "C" {
     const int jcols = nExog+nEndo+nsPred+nsForw; // Num of Jacobian columns
 
     mxFldp = mxGetField(M_, 0, "var_order_endo_names");
-    const int nendo = (int) mxGetM(mxFldp);
-    const int widthEndo = (int) mxGetN(mxFldp);
+    const int nendo = static_cast<int>(mxGetM(mxFldp));
+    const int widthEndo = static_cast<int>(mxGetN(mxFldp));
     std::vector<std::string> endoNames;
     DynareMxArrayToString(mxFldp, nendo, widthEndo, endoNames);
 
     mxFldp = mxGetField(M_, 0, "exo_names");
-    const int nexo = (int) mxGetM(mxFldp);
-    const int widthExog = (int) mxGetN(mxFldp);
+    const int nexo = static_cast<int>(mxGetM(mxFldp));
+    const int widthExog = static_cast<int>(mxGetN(mxFldp));
     std::vector<std::string> exoNames;
     DynareMxArrayToString(mxFldp, nexo, widthExog, exoNames);
 
-    if ((nEndo != nendo) || (nExog != nexo))
+    if (nEndo != nendo || nExog != nexo)
       DYN_MEX_FUNC_ERR_MSG_TXT("Incorrect number of input parameters.");
 
-    TwoDMatrix *g1m = nullptr;
-    TwoDMatrix *g2m = nullptr;
-    TwoDMatrix *g3m = nullptr;
+    std::unique_ptr<TwoDMatrix> g1m, g2m, g3m;
     // derivatives passed as arguments */
     if (nrhs > 3)
       {
         const mxArray *g1 = prhs[3];
-        int m = (int) mxGetM(g1);
-        int n = (int) mxGetN(g1);
-        g1m = new TwoDMatrix(m, n, Vector{ConstVector{g1}});
+        int m = static_cast<int>(mxGetM(g1));
+        int n = static_cast<int>(mxGetN(g1));
+        g1m = std::make_unique<TwoDMatrix>(m, n, Vector{ConstVector{g1}});
         if (nrhs > 4)
           {
             const mxArray *g2 = prhs[4];
-            int m = (int) mxGetM(g2);
-            int n = (int) mxGetN(g2);
-            g2m = new TwoDMatrix(m, n, Vector{ConstVector{g2}});
+            int m = static_cast<int>(mxGetM(g2));
+            int n = static_cast<int>(mxGetN(g2));
+            g2m = std::make_unique<TwoDMatrix>(m, n, Vector{ConstVector{g2}});
             if (nrhs > 5)
               {
                 const mxArray *g3 = prhs[5];
-                int m = (int) mxGetM(g3);
-                int n = (int) mxGetN(g3);
-                g3m = new TwoDMatrix(m, n, Vector{ConstVector{g3}});
+                int m = static_cast<int>(mxGetM(g3));
+                int n = static_cast<int>(mxGetN(g3));
+                g3m = std::make_unique<TwoDMatrix>(m, n, Vector{ConstVector{g3}});
               }
           }
       }
@@ -222,13 +210,10 @@ extern "C" {
 
     try
       {
-        // make journal name and journal
-        std::string jName(fName); //params.basename);
-        jName += ".jnl";
-        Journal journal(jName.c_str());
+        Journal journal(fName + ".jnl");
 
         std::unique_ptr<DynamicModelAC> dynamicModelFile;
-        if (use_dll == 1)
+        if (use_dll)
           dynamicModelFile = std::make_unique<DynamicModelDLL>(fName);
         else
           dynamicModelFile = std::make_unique<DynamicModelMFile>(fName);
@@ -241,17 +226,17 @@ extern "C" {
                            ySteady, vCov, modParams, nStat, nPred, nForw, nBoth,
                            jcols, NNZD, nSteps, kOrder, journal, std::move(dynamicModelFile),
                            sstol, var_order_vp, llincidence, qz_criterium,
-                           g1m, g2m, g3m);
+                           std::move(g1m), std::move(g2m), std::move(g3m));
 
         // construct main K-order approximation class
 
-        Approximation app(dynare, journal,  nSteps, false, qz_criterium);
+        Approximation app(dynare, journal, nSteps, false, qz_criterium);
         // run stochastic steady
         app.walkStochSteady();
 
         /* Write derivative outputs into memory map */
         std::map<std::string, ConstTwoDMatrix> mm;
-        app.getFoldDecisionRule().writeMMap(mm, std::string());
+        app.getFoldDecisionRule().writeMMap(mm, "");
 
         // get latest ysteady
         ySteady = dynare.getSteady();
@@ -261,25 +246,24 @@ extern "C" {
             /* Set the output pointer to the output matrix ysteady. */
             auto cit = mm.begin();
             ++cit;
-            plhs[1] = mxCreateDoubleMatrix((*cit).second.numRows(), (*cit).second.numCols(), mxREAL);
+            plhs[1] = mxCreateDoubleMatrix(cit->second.numRows(), cit->second.numCols(), mxREAL);
 
             // Copy Dynare++ matrix into MATLAB matrix
-            const ConstVector &vec = (*cit).second.getData();
+            const ConstVector &vec = cit->second.getData();
             assert(vec.skip() == 1);
-            memcpy(mxGetPr(plhs[1]), vec.base(), vec.length() * sizeof(double));
+            std::copy_n(vec.base(), vec.length(), mxGetPr(plhs[1]));
           }
         if (kOrder >= 2)
           {
             int ii = 1;
-            for (auto cit = mm.begin();
-                 ((cit != mm.end()) && (ii < nlhs)); ++cit)
+            for (auto cit = mm.begin(); cit != mm.end() && ii < nlhs; ++cit)
               {
-                plhs[ii] = mxCreateDoubleMatrix((*cit).second.numRows(), (*cit).second.numCols(), mxREAL);
+                plhs[ii] = mxCreateDoubleMatrix(cit->second.numRows(), cit->second.numCols(), mxREAL);
 
                 // Copy Dynare++ matrix into MATLAB matrix
-                const ConstVector &vec = (*cit).second.getData();
+                const ConstVector &vec = cit->second.getData();
                 assert(vec.skip() == 1);
-                memcpy(mxGetPr(plhs[ii]), vec.base(), vec.length() * sizeof(double));
+                std::copy_n(vec.base(), vec.length(), mxGetPr(plhs[ii]));
 
                 ++ii;
 
@@ -312,9 +296,7 @@ extern "C" {
     catch (const KordException &e)
       {
         e.print();
-        std::ostringstream strstrm;
-        strstrm << "dynare:k_order_perturbation: Caught Kord exception: " << e.get_message();
-        DYN_MEX_FUNC_ERR_MSG_TXT(strstrm.str().c_str());
+        DYN_MEX_FUNC_ERR_MSG_TXT("dynare:k_order_perturbation: Caught Kord exception: " + e.get_message());
       }
     catch (const TLException &e)
       {
@@ -328,17 +310,13 @@ extern "C" {
       }
     catch (const DynareException &e)
       {
-        std::ostringstream strstrm;
-        strstrm << "dynare:k_order_perturbation: Caught KordDynare exception: " << e.message();
-        DYN_MEX_FUNC_ERR_MSG_TXT(strstrm.str().c_str());
+        DYN_MEX_FUNC_ERR_MSG_TXT("dynare:k_order_perturbation: Caught KordDynare exception: " + e.message());
       }
     catch (const ogu::Exception &e)
       {
-        std::ostringstream strstrm;
-        strstrm << "dynare:k_order_perturbation: Caught general exception: " << e.message();
-        DYN_MEX_FUNC_ERR_MSG_TXT(strstrm.str().c_str());
+        DYN_MEX_FUNC_ERR_MSG_TXT("dynare:k_order_perturbation: Caught general exception: " + e.message());
       }
     plhs[0] = mxCreateDoubleScalar(0);
   } // end of mexFunction()
 } // end of extern C
-#endif // ifdef MATLAB_MEX_FILE  to exclude mexFunction for other applications
+