diff --git a/mex/build/estimation.am b/mex/build/estimation.am
index 3d919c7de..b2a92a326 100644
--- a/mex/build/estimation.am
+++ b/mex/build/estimation.am
@@ -63,6 +63,5 @@ nodist_logMHMCMCposterior_SOURCES = \
 	$(COMMON_SRCS) \
 	$(TOPDIR)/Proposal.cc \
 	$(TOPDIR)/Proposal.hh \
-	$(TOPDIR)/RandomWalkMetropolisHastings.cc \
 	$(TOPDIR)/RandomWalkMetropolisHastings.hh \
 	$(TOPDIR)/logMHMCMCposterior.cc
diff --git a/mex/sources/estimation/InitializeKalmanFilter.cc b/mex/sources/estimation/InitializeKalmanFilter.cc
index c409b23d0..182b6c24b 100644
--- a/mex/sources/estimation/InitializeKalmanFilter.cc
+++ b/mex/sources/estimation/InitializeKalmanFilter.cc
@@ -55,31 +55,6 @@ InitializeKalmanFilter::InitializeKalmanFilter(const std::string &dynamicDllFile
                              zeta_back_mixed[i]) - zeta_varobs_back_mixed.begin());
 
 }
-// initialise parameter dependent KF matrices only but not Ps
-void
-InitializeKalmanFilter::initialize(VectorView &steadyState, const Vector &deepParams, Matrix &R,
-                                   const Matrix &Q, Matrix &RQRt, Matrix &T,
-                                   double &penalty, const MatrixConstView &dataView,
-                                   MatrixView &detrendedDataView, int &info)
-{
-  modelSolution.compute(steadyState, deepParams, g_x, g_u);
-  detrendData.detrend(steadyState, dataView, detrendedDataView);
-
-  setT(T, info);
-  setRQR(R, Q, RQRt, info);
-}
-
-// initialise all KF matrices
-void
-InitializeKalmanFilter::initialize(VectorView &steadyState, const Vector &deepParams, Matrix &R,
-                                   const Matrix &Q, Matrix &RQRt, Matrix &T, Matrix &Pstar, Matrix &Pinf,
-                                   double &penalty, const MatrixConstView &dataView,
-                                   MatrixView &detrendedDataView, int &info)
-{
-  initialize(steadyState, deepParams, R, Q, RQRt, T, penalty, dataView, detrendedDataView, info);
-  setPstar(Pstar, Pinf, T, RQRt, info);
-}
-
 void
 InitializeKalmanFilter::setT(Matrix &T, int &info)
 {
@@ -89,7 +64,7 @@ InitializeKalmanFilter::setT(Matrix &T, int &info)
 }
 
 void
-InitializeKalmanFilter::setRQR(Matrix &R, const Matrix &Q, Matrix &RQRt, int &info)
+InitializeKalmanFilter::setRQR(Matrix &R, const MatrixView &Q, Matrix &RQRt, int &info)
 {
   mat::assignByVectors(R, mat::nullVec, mat::nullVec, g_u, zeta_varobs_back_mixed, mat::nullVec);
 
diff --git a/mex/sources/estimation/InitializeKalmanFilter.hh b/mex/sources/estimation/InitializeKalmanFilter.hh
index b8d7f6cba..399934255 100644
--- a/mex/sources/estimation/InitializeKalmanFilter.hh
+++ b/mex/sources/estimation/InitializeKalmanFilter.hh
@@ -48,12 +48,30 @@ public:
                          const std::vector<size_t> &zeta_varobs_back_mixed_arg,
                          double qz_criterium_arg, double lyapunov_tol_arg, int &info);
   virtual ~InitializeKalmanFilter();
-  // initialise all KF matrices
-  void initialize(VectorView &steadyState, const Vector &deepParams, Matrix &R, const Matrix &Q, Matrix &RQRt,
-                  Matrix &T, Matrix &Pstar, Matrix &Pinf, double &penalty, const MatrixConstView &dataView, MatrixView &detrendedDataView, int &info);
   // initialise parameter dependent KF matrices only but not Ps
-  void initialize(VectorView &steadyState, const Vector &deepParams, Matrix &R, const Matrix &Q, Matrix &RQRt,
-                  Matrix &T, double &penalty, const MatrixConstView &dataView, MatrixView &detrendedDataView, int &info);
+  template <class VEC>
+  void initialize(VEC &steadyState, const VectorView &deepParams, Matrix &R,
+				     const MatrixView &Q, Matrix &RQRt, Matrix &T,
+				     double &penalty, const MatrixConstView &dataView,
+				     MatrixView &detrendedDataView, int &info)
+  {
+    modelSolution.compute(steadyState, deepParams, g_x, g_u);
+    detrendData.detrend(steadyState, dataView, detrendedDataView);
+
+    setT(T, info);
+    setRQR(R, Q, RQRt, info);
+  }
+
+  // initialise all KF matrices
+  template <class VEC>
+  void initialize(VEC &steadyState, const VectorView &deepParams, Matrix &R,
+				     const MatrixView &Q, Matrix &RQRt, Matrix &T, Matrix &Pstar, Matrix &Pinf,
+				     double &penalty, const MatrixConstView &dataView,
+				     MatrixView &detrendedDataView, int &info)
+  {
+    initialize(steadyState, deepParams, R, Q, RQRt, T, penalty, dataView, detrendedDataView, info);
+    setPstar(Pstar, Pinf, T, RQRt, info);
+  }
 
 private:
   const double lyapunov_tol;
@@ -68,7 +86,7 @@ private:
   Matrix g_u;
   Matrix Rt, RQ;
   void setT(Matrix &T, int &info);
-  void setRQR(Matrix &R, const Matrix &Q, Matrix &RQRt, int &info);
+  void setRQR(Matrix &R, const MatrixView &Q, Matrix &RQRt, int &info);
   void setPstar(Matrix &Pstar, Matrix &Pinf, const Matrix &T, const Matrix &RQRt, int &info);
 
 };
diff --git a/mex/sources/estimation/KalmanFilter.cc b/mex/sources/estimation/KalmanFilter.cc
index 914f7b542..77e8d0804 100644
--- a/mex/sources/estimation/KalmanFilter.cc
+++ b/mex/sources/estimation/KalmanFilter.cc
@@ -73,36 +73,6 @@ KalmanFilter::compute_zeta_varobs_back_mixed(const std::vector<size_t> &zeta_bac
   return zeta_varobs_back_mixed;
 }
 
-double
-KalmanFilter::compute(const MatrixConstView &dataView, VectorView &steadyState,
-                      const Matrix &Q, const Matrix &H, const Vector &deepParams,
-                      VectorView &vll, MatrixView &detrendedDataView,
-                      size_t start, size_t period, double &penalty, int &info)
-{
-  double lik = INFINITY;
-  try
-    {
-      if (period == 0) // initialise all KF matrices
-        initKalmanFilter.initialize(steadyState, deepParams, R, Q, RQRt, T, Pstar, Pinf,
-                                    penalty, dataView, detrendedDataView, info);
-      else             // initialise parameter dependent KF matrices only but not Ps
-        initKalmanFilter.initialize(steadyState, deepParams, R, Q, RQRt, T,
-                                    penalty, dataView, detrendedDataView, info);
-
-      lik = filter(detrendedDataView, H, vll, start, info);
-    }
-  catch (const DecisionRules::BlanchardKahnException &bke)
-    {
-      info = 22;
-      return penalty;
-    }
-
-  if (info != 0)
-    return penalty;
-  else
-    return lik;
-
-};
 
 /**
  * Multi-variate standard Kalman Filter
@@ -110,7 +80,7 @@ KalmanFilter::compute(const MatrixConstView &dataView, VectorView &steadyState,
 double
 KalmanFilter::filter(const MatrixView &detrendedDataView,  const Matrix &H, VectorView &vll, size_t start, int &info)
 {
-  double loglik = 0.0, ll, logFdet, Fdet, dvtFinvVt;
+  double loglik = 0.0, ll, logFdet = 0.0, Fdet, dvtFinvVt;
   size_t p = Finv.getRows();
   bool nonstationary = true;
   a_init.setAll(0.0);
diff --git a/mex/sources/estimation/KalmanFilter.hh b/mex/sources/estimation/KalmanFilter.hh
index 89b8b9d61..64fb03a86 100644
--- a/mex/sources/estimation/KalmanFilter.hh
+++ b/mex/sources/estimation/KalmanFilter.hh
@@ -54,10 +54,36 @@ public:
                double qz_criterium_arg, const std::vector<size_t> &varobs_arg,
                double riccati_tol_arg, double lyapunov_tol_arg, int &info);
 
-  double compute(const MatrixConstView &dataView, VectorView &steadyState,
-                 const Matrix &Q, const Matrix &H, const Vector &deepParams,
+  template <class VEC>
+  double compute(const MatrixConstView &dataView, VEC &steadyState,
+                 const MatrixView &Q, const Matrix &H, const VectorView &deepParams,
                  VectorView &vll, MatrixView &detrendedDataView, size_t start, size_t period,
-                 double &penalty, int &info);
+                 double &penalty, int &info)
+  {
+    double lik = INFINITY;
+    try
+      {
+	if (period == 0) // initialise all KF matrices
+	  initKalmanFilter.initialize(steadyState, deepParams, R, Q, RQRt, T, Pstar, Pinf,
+				      penalty, dataView, detrendedDataView, info);
+	else             // initialise parameter dependent KF matrices only but not Ps
+	  initKalmanFilter.initialize(steadyState, deepParams, R, Q, RQRt, T,
+				      penalty, dataView, detrendedDataView, info);
+
+	lik = filter(detrendedDataView, H, vll, start, info);
+      }
+    catch (const DecisionRules::BlanchardKahnException &bke)
+      {
+	info = 22;
+	return penalty;
+      }
+
+    if (info != 0)
+      return penalty;
+    else
+      return lik;
+
+  };
 
 private:
   const std::vector<size_t> zeta_varobs_back_mixed;
diff --git a/mex/sources/estimation/LogLikelihoodMain.cc b/mex/sources/estimation/LogLikelihoodMain.cc
index 5edcb3513..5371acd04 100644
--- a/mex/sources/estimation/LogLikelihoodMain.cc
+++ b/mex/sources/estimation/LogLikelihoodMain.cc
@@ -44,23 +44,4 @@ LogLikelihoodMain::~LogLikelihoodMain()
 
 }
 
-double
-LogLikelihoodMain::compute(Matrix &steadyState, const Vector &estParams, Vector &deepParams, const MatrixConstView &data, Matrix &Q, Matrix &H, size_t start, int &info)
-{
-  double logLikelihood = 0;
-  for (size_t i = 0; i < estSubsamples.size(); ++i)
-    {
-      VectorView vSteadyState = mat::get_col(steadyState, i);
-
-      MatrixConstView dataView(data, 0, estSubsamples[i].startPeriod,
-                               data.getRows(), estSubsamples[i].endPeriod-estSubsamples[i].startPeriod+1);
-      MatrixView detrendedDataView(detrendedData, 0, estSubsamples[i].startPeriod,
-                                   data.getRows(), estSubsamples[i].endPeriod-estSubsamples[i].startPeriod+1);
-
-      VectorView vllView(vll, estSubsamples[i].startPeriod, estSubsamples[i].endPeriod-estSubsamples[i].startPeriod+1);
-      logLikelihood += logLikelihoodSubSample.compute(vSteadyState, dataView, estParams, deepParams,
-                                                      Q, H, vllView, detrendedDataView, info, start, i);
-    }
-  return logLikelihood;
-};
 
diff --git a/mex/sources/estimation/LogLikelihoodMain.hh b/mex/sources/estimation/LogLikelihoodMain.hh
index 1ab2c8fb8..546e122df 100644
--- a/mex/sources/estimation/LogLikelihoodMain.hh
+++ b/mex/sources/estimation/LogLikelihoodMain.hh
@@ -50,14 +50,29 @@ public:
    * Q and H KF matrices of shock and measurement error varinaces and covariances
    * KF logLikelihood calculation start period.
    */
-  double compute(Matrix &steadyState, const Vector &estParams, Vector &deepParams, const MatrixConstView &data,
-                 Matrix &Q, Matrix &H, size_t presampleStart, int &info); // for calls from estimation and to set Steady State
 
-  Vector &
-  getVll()
+  template <class VEC1, class VEC2>
+  double compute(VEC1 &steadyState, VEC2 &estParams, VectorView &deepParams, const MatrixConstView &data, 
+		 MatrixView &Q, Matrix &H, size_t start, int &info)
   {
-    return vll;
+    double logLikelihood = 0;
+    for (size_t i = 0; i < estSubsamples.size(); ++i)
+      {
+	VectorView vSteadyState (steadyState,0,steadyState.getSize());
+
+	MatrixConstView dataView(data, 0, estSubsamples[i].startPeriod,
+				 data.getRows(), estSubsamples[i].endPeriod-estSubsamples[i].startPeriod+1);
+	MatrixView detrendedDataView(detrendedData, 0, estSubsamples[i].startPeriod,
+				     data.getRows(), estSubsamples[i].endPeriod-estSubsamples[i].startPeriod+1);
+
+	VectorView vllView(vll, estSubsamples[i].startPeriod, estSubsamples[i].endPeriod-estSubsamples[i].startPeriod+1);
+	logLikelihood += logLikelihoodSubSample.compute(vSteadyState, dataView, estParams, deepParams,
+							Q, H, vllView, detrendedDataView, info, start, i);
+      }
+    return logLikelihood;
   };
+
+  Vector &getVll() { return vll; };
 };
 
 #endif // !defined(E126AEF5_AC28_400a_821A_3BCFD1BC4C22__INCLUDED_)
diff --git a/mex/sources/estimation/LogLikelihoodSubSample.cc b/mex/sources/estimation/LogLikelihoodSubSample.cc
index c601a59c6..4c5063db2 100644
--- a/mex/sources/estimation/LogLikelihoodSubSample.cc
+++ b/mex/sources/estimation/LogLikelihoodSubSample.cc
@@ -25,8 +25,6 @@
 
 //#include "LogLikelihoodSubSample.hh"
 #include "LogLikelihoodMain.hh" // use ...Main.hh for testing only
-#include <algorithm>
-#include "LapackBindings.hh"
 
 LogLikelihoodSubSample::~LogLikelihoodSubSample()
 {
@@ -42,133 +40,3 @@ LogLikelihoodSubSample::LogLikelihoodSubSample(const std::string &dynamicDllFile
 {
 };
 
-double
-LogLikelihoodSubSample::compute(VectorView &steadyState, const MatrixConstView &dataView, const Vector &estParams, Vector &deepParams,
-                                Matrix &Q, Matrix &H, VectorView &vll, MatrixView &detrendedDataView, int &info, size_t start, size_t period)
-{
-  penalty = startPenalty;
-  logLikelihood = startPenalty;
-
-  updateParams(estParams, deepParams, Q, H, period);
-  if (info == 0)
-    logLikelihood = kalmanFilter.compute(dataView, steadyState,  Q, H, deepParams, vll, detrendedDataView, start, period, penalty,  info);
-  //  else
-  //    logLikelihood+=penalty;
-
-  return logLikelihood;
-
-};
-
-void
-LogLikelihoodSubSample::updateParams(const Vector &estParams, Vector &deepParams,
-                                     Matrix &Q, Matrix &H, size_t period)
-{
-  size_t i, k, k1, k2;
-  int test;
-  bool found;
-  std::vector<size_t>::const_iterator it;
-  info = 0;
-
-  for (i = 0; i <  estParams.getSize(); ++i)
-    {
-      found = false;
-      it = find(estiParDesc.estParams[i].subSampleIDs.begin(),
-                estiParDesc.estParams[i].subSampleIDs.end(), period);
-      if (it != estiParDesc.estParams[i].subSampleIDs.end())
-        found = true;
-      if (found)
-        {
-          switch (estiParDesc.estParams[i].ptype)
-            {
-            case EstimatedParameter::shock_SD:
-              k = estiParDesc.estParams[i].ID1;
-              Q(k, k) = estParams(i)*estParams(i);
-              break;
-
-            case EstimatedParameter::measureErr_SD:
-              k = estiParDesc.estParams[i].ID1;
-              H(k, k) = estParams(i)*estParams(i);
-              break;
-
-            case EstimatedParameter::shock_Corr:
-              k1 = estiParDesc.estParams[i].ID1;
-              k2 = estiParDesc.estParams[i].ID2;
-              Q(k1, k2) = estParams(i)*sqrt(Q(k1, k1)*Q(k2, k2));
-              Q(k2, k1) = Q(k1, k2);
-              //   [CholQ,testQ] = chol(Q);
-              test = lapack::choleskyDecomp(Q, "L");
-              if (test > 0)
-                {
-                  mexPrintf("Caugth unhandled exception with cholesky of Q matrix: ");
-                  logLikelihood = penalty;
-                  info = 1;
-                }
-              else if (test < 0)
-                {
-                  // The variance-covariance matrix of the structural innovations is not definite positive.
-                  // We have to compute the eigenvalues of this matrix in order to build the penalty.
-                  double delta = 0;
-                  eigQ.calculate(Q);  // get eigenvalues
-                  //k = find(a < 0);
-                  if (eigQ.hasConverged())
-                    {
-                      const Vector &evQ = eigQ.getD();
-                      for (i = 0; i < evQ.getSize(); ++i)
-                        if (evQ(i) < 0)
-                          delta -= evQ(i);
-                    }
-
-                  logLikelihood = penalty+delta;
-                  info = 43;
-                } // if
-              break;
-
-            case EstimatedParameter::measureErr_Corr:
-              k1 = estiParDesc.estParams[i].ID1;
-              k2 = estiParDesc.estParams[i].ID2;
-              //      H(k1,k2) = xparam1(i)*sqrt(H(k1,k1)*H(k2,k2));
-              //      H(k2,k1) = H(k1,k2);
-              H(k1, k2) = estParams(i)*sqrt(H(k1, k1)*H(k2, k2));
-              H(k2, k1) = H(k1, k2);
-
-              //[CholH,testH] = chol(H);
-              test = lapack::choleskyDecomp(H, "L");
-              if (test > 0)
-                {
-                  mexPrintf("Caugth unhandled exception with cholesky of Q matrix: ");
-                  logLikelihood = penalty;
-                  info = 1;
-                }
-              else if (test < 0)
-                {
-                  // The variance-covariance matrix of the measurement errors is not definite positive.
-                  // We have to compute the eigenvalues of this matrix in order to build the penalty.
-                  //a = diag(eig(H));
-                  double delta = 0;
-                  eigH.calculate(H);  // get eigenvalues
-                  //k = find(a < 0);
-                  if (eigH.hasConverged())
-                    {
-                      const Vector &evH = eigH.getD();
-                      for (i = 0; i < evH.getSize(); ++i)
-                        if (evH(i) < 0)
-                          delta -= evH(i);
-                    }
-                  logLikelihood = penalty+delta;
-                  info = 44;
-                } //   end if
-              break;
-
-              //if estim_params_.np > 0  // i.e. num of deep parameters >0
-            case EstimatedParameter::deepPar:
-              k = estiParDesc.estParams[i].ID1;
-              deepParams(k) = estParams(i);
-              break;
-            default:
-              logLikelihood = penalty;
-              info = 1;
-            } // end switch
-        } // end found
-    } //end for
-};
-
diff --git a/mex/sources/estimation/LogLikelihoodSubSample.hh b/mex/sources/estimation/LogLikelihoodSubSample.hh
index dec031af2..9df4f5fd6 100644
--- a/mex/sources/estimation/LogLikelihoodSubSample.hh
+++ b/mex/sources/estimation/LogLikelihoodSubSample.hh
@@ -26,9 +26,11 @@
 #if !defined(DF8B7AF5_8169_4587_9037_2CD2C82E2DDF__INCLUDED_)
 #define DF8B7AF5_8169_4587_9037_2CD2C82E2DDF__INCLUDED_
 
+#include <algorithm>
 #include "EstimatedParametersDescription.hh"
 #include "KalmanFilter.hh"
 #include "VDVEigDecomposition.hh"
+#include "LapackBindings.hh"
 
 class LogLikelihoodSubSample
 {
@@ -39,8 +41,23 @@ public:
                          const std::vector<size_t> &zeta_mixed_arg, const std::vector<size_t> &zeta_static_arg, const double qz_criterium,
                          const std::vector<size_t> &varobs_arg, double riccati_tol_in, double lyapunov_tol, int &info);
 
-  double compute(VectorView &steadyState, const MatrixConstView &dataView, const Vector &estParams, Vector &deepParams,
-                 Matrix &Q, Matrix &H, VectorView &vll, MatrixView &detrendedDataView, int &info,  size_t start, size_t period);
+  template <class VEC1, class VEC2>
+  double compute(VEC1 &steadyState, const MatrixConstView &dataView, VEC2 &estParams, VectorView &deepParams,
+		 MatrixView &Q, Matrix &H, VectorView &vll, MatrixView &detrendedDataView, int &info, size_t start, size_t period)
+  {
+    penalty = startPenalty;
+    logLikelihood = startPenalty;
+
+    updateParams(estParams, deepParams, Q, H, period);
+    if (info == 0)
+      logLikelihood = kalmanFilter.compute(dataView, steadyState,  Q, H, deepParams, vll, detrendedDataView, start, period, penalty,  info);
+    //  else
+    //    logLikelihood+=penalty;
+
+    return logLikelihood;
+
+  };
+
   virtual ~LogLikelihoodSubSample();
 
 private:
@@ -53,8 +70,119 @@ private:
   int &info;
 
   // methods
-  void updateParams(const Vector &estParams, Vector &deepParams,
-                    Matrix &Q, Matrix &H, size_t period);
+  template <class VEC>
+  void updateParams(VEC &estParams, VectorView &deepParams,
+                    MatrixView &Q, Matrix &H, size_t period)
+  {
+    size_t i, k, k1, k2;
+    int test;
+    bool found;
+    std::vector<size_t>::const_iterator it;
+    info = 0;
+
+    for (i = 0; i <  estParams.getSize(); ++i)
+      {
+	found = false;
+	it = find(estiParDesc.estParams[i].subSampleIDs.begin(),
+		  estiParDesc.estParams[i].subSampleIDs.end(), period);
+	if (it != estiParDesc.estParams[i].subSampleIDs.end())
+	  found = true;
+	if (found)
+	  {
+	    switch (estiParDesc.estParams[i].ptype)
+	      {
+	      case EstimatedParameter::shock_SD:
+		k = estiParDesc.estParams[i].ID1;
+		Q(k, k) = estParams(i)*estParams(i);
+		break;
+
+	      case EstimatedParameter::measureErr_SD:
+		k = estiParDesc.estParams[i].ID1;
+		H(k, k) = estParams(i)*estParams(i);
+		break;
+
+	      case EstimatedParameter::shock_Corr:
+		k1 = estiParDesc.estParams[i].ID1;
+		k2 = estiParDesc.estParams[i].ID2;
+		Q(k1, k2) = estParams(i)*sqrt(Q(k1, k1)*Q(k2, k2));
+		Q(k2, k1) = Q(k1, k2);
+		//   [CholQ,testQ] = chol(Q);
+		test = lapack::choleskyDecomp(Q, "L");
+		if (test > 0)
+		  {
+		    mexPrintf("Caugth unhandled exception with cholesky of Q matrix: ");
+		    logLikelihood = penalty;
+		    info = 1;
+		  }
+		else if (test < 0)
+		  {
+		    // The variance-covariance matrix of the structural innovations is not definite positive.
+		    // We have to compute the eigenvalues of this matrix in order to build the penalty.
+		    double delta = 0;
+		    eigQ.calculate(Q);  // get eigenvalues
+		    //k = find(a < 0);
+		    if (eigQ.hasConverged())
+		      {
+			const Vector &evQ = eigQ.getD();
+			for (i = 0; i < evQ.getSize(); ++i)
+			  if (evQ(i) < 0)
+			    delta -= evQ(i);
+		      }
+
+		    logLikelihood = penalty+delta;
+		    info = 43;
+		  } // if
+		break;
+
+	      case EstimatedParameter::measureErr_Corr:
+		k1 = estiParDesc.estParams[i].ID1;
+		k2 = estiParDesc.estParams[i].ID2;
+		//      H(k1,k2) = xparam1(i)*sqrt(H(k1,k1)*H(k2,k2));
+		//      H(k2,k1) = H(k1,k2);
+		H(k1, k2) = estParams(i)*sqrt(H(k1, k1)*H(k2, k2));
+		H(k2, k1) = H(k1, k2);
+
+		//[CholH,testH] = chol(H);
+		test = lapack::choleskyDecomp(H, "L");
+		if (test > 0)
+		  {
+		    mexPrintf("Caugth unhandled exception with cholesky of Q matrix: ");
+		    logLikelihood = penalty;
+		    info = 1;
+		  }
+		else if (test < 0)
+		  {
+		    // The variance-covariance matrix of the measurement errors is not definite positive.
+		    // We have to compute the eigenvalues of this matrix in order to build the penalty.
+		    //a = diag(eig(H));
+		    double delta = 0;
+		    eigH.calculate(H);  // get eigenvalues
+		    //k = find(a < 0);
+		    if (eigH.hasConverged())
+		      {
+			const Vector &evH = eigH.getD();
+			for (i = 0; i < evH.getSize(); ++i)
+			  if (evH(i) < 0)
+			    delta -= evH(i);
+		      }
+		    logLikelihood = penalty+delta;
+		    info = 44;
+		  } //   end if
+		break;
+
+		//if estim_params_.np > 0  // i.e. num of deep parameters >0
+	      case EstimatedParameter::deepPar:
+		k = estiParDesc.estParams[i].ID1;
+		deepParams(k) = estParams(i);
+		break;
+	      default:
+		logLikelihood = penalty;
+		info = 1;
+	      } // end switch
+	  } // end found
+      } //end for
+  };
+
 
 };
 
diff --git a/mex/sources/estimation/LogPosteriorDensity.cc b/mex/sources/estimation/LogPosteriorDensity.cc
index f4436411e..b99e0b8ee 100644
--- a/mex/sources/estimation/LogPosteriorDensity.cc
+++ b/mex/sources/estimation/LogPosteriorDensity.cc
@@ -40,12 +40,6 @@ LogPosteriorDensity::LogPosteriorDensity(const std::string &modName, EstimatedPa
 
 }
 
-double
-LogPosteriorDensity::compute(Matrix &steadyState, const Vector &estParams, Vector &deepParams, const MatrixConstView &data, Matrix &Q, Matrix &H, size_t presampleStart, int &info)
-{
-  return -logLikelihoodMain.compute(steadyState, estParams, deepParams, data, Q, H, presampleStart, info)
-    -logPriorDensity.compute(estParams);
-}
 
 /**
  * vector of log likelihoods for each Kalman step
diff --git a/mex/sources/estimation/LogPosteriorDensity.hh b/mex/sources/estimation/LogPosteriorDensity.hh
index 53e2aa401..5c378c822 100644
--- a/mex/sources/estimation/LogPosteriorDensity.hh
+++ b/mex/sources/estimation/LogPosteriorDensity.hh
@@ -48,8 +48,15 @@ public:
                       const std::vector<size_t> &zeta_fwrd_arg, const std::vector<size_t> &zeta_back_arg, const std::vector<size_t> &zeta_mixed_arg,
                       const std::vector<size_t> &zeta_static_arg, const double qz_criterium_arg, const std::vector<size_t> &varobs_arg,
                       double riccati_tol_arg, double lyapunov_tol_arg, int &info_arg);
-  double compute(Matrix &steadyState, const Vector &estParams, Vector &deepParams,
-                 const MatrixConstView &data, Matrix &Q, Matrix &H, size_t presampleStart, int &info);
+
+  template <class VEC1, class VEC2>
+  double
+  compute(VEC1 &steadyState, VEC2 &estParams, VectorView &deepParams, const MatrixConstView &data, MatrixView &Q, Matrix &H, size_t presampleStart, int &info)
+  {
+    return -logLikelihoodMain.compute(steadyState, estParams, deepParams, data, Q, H, presampleStart, info)
+      -logPriorDensity.compute(estParams);
+  }
+
   Vector&getLikVector();
 
 };
diff --git a/mex/sources/estimation/LogPriorDensity.cc b/mex/sources/estimation/LogPriorDensity.cc
index e23eeb751..55364d5c7 100644
--- a/mex/sources/estimation/LogPriorDensity.cc
+++ b/mex/sources/estimation/LogPriorDensity.cc
@@ -33,20 +33,6 @@ LogPriorDensity::LogPriorDensity(EstimatedParametersDescription &estParsDesc_arg
 {
 };
 
-double
-LogPriorDensity::compute(const Vector &ep)
-{
-  assert(estParsDesc.estParams.size() == ep.getSize());
-  double logPriorDensity = 0;
-  for (size_t i = 0; i <  ep.getSize(); ++i)
-    {
-      logPriorDensity += log(((*(estParsDesc.estParams[i]).prior)).pdf(ep(i)));
-      if (std::isinf(fabs(logPriorDensity)))
-        return logPriorDensity;
-    }
-  return logPriorDensity;
-};
-
 /**
  * Return random number for prior fromits distribution
  */
diff --git a/mex/sources/estimation/LogPriorDensity.hh b/mex/sources/estimation/LogPriorDensity.hh
index d74fa9825..ef0858e52 100644
--- a/mex/sources/estimation/LogPriorDensity.hh
+++ b/mex/sources/estimation/LogPriorDensity.hh
@@ -37,7 +37,20 @@ public:
   LogPriorDensity(EstimatedParametersDescription &estParsDesc);
   virtual ~LogPriorDensity();
 
-  double compute(const Vector &estParams);
+  template<class VEC>
+  double compute(VEC &ep)
+  {
+    assert(estParsDesc.estParams.size() == ep.getSize());
+    double logPriorDensity = 0;
+    for (size_t i = 0; i <  ep.getSize(); ++i)
+      {
+	logPriorDensity += log(((*(estParsDesc.estParams[i]).prior)).pdf(ep(i)));
+	if (std::isinf(fabs(logPriorDensity)))
+	  return logPriorDensity;
+      }
+    return logPriorDensity;
+  };
+
   void computeNewParams(Vector &newParams);
 
 private:
diff --git a/mex/sources/estimation/ModelSolution.cc b/mex/sources/estimation/ModelSolution.cc
index 3ba3f8268..4eb8842a4 100644
--- a/mex/sources/estimation/ModelSolution.cc
+++ b/mex/sources/estimation/ModelSolution.cc
@@ -51,41 +51,4 @@ ModelSolution::ModelSolution(const std::string &dynamicDllFile,  size_t n_endo_a
             back_inserter(zeta_back_mixed));
 }
 
-void
-ModelSolution::compute(VectorView &steadyState, const Vector &deepParams, Matrix &ghx, Matrix &ghu) throw (DecisionRules::BlanchardKahnException, GeneralizedSchurDecomposition::GSDException)
-{
-  // compute Steady State
-  ComputeSteadyState(steadyState, deepParams);
-
-  // then get jacobian and
-
-  ComputeModelSolution(steadyState, deepParams, ghx, ghu);
-
-}
-
-void
-ModelSolution::ComputeModelSolution(VectorView &steadyState, const Vector &deepParams, Matrix &ghx, Matrix &ghu) throw (DecisionRules::BlanchardKahnException, GeneralizedSchurDecomposition::GSDException)
-{
-  // set extended Steady State
-
-  for (size_t i = 0; i < zeta_back_mixed.size(); i++)
-    llXsteadyState(i) = steadyState(zeta_back_mixed[i]);
-
-  for (size_t i = 0; i < n_endo; i++)
-    llXsteadyState(zeta_back_mixed.size() + i) = steadyState(i);
-
-  for (size_t i = 0; i < zeta_fwrd_mixed.size(); i++)
-    llXsteadyState(zeta_back_mixed.size() + n_endo + i) = steadyState(zeta_fwrd_mixed[i]);
-
-  //get jacobian
-  dynamicDLLp.eval(llXsteadyState, Mx, deepParams, steadyState, residual, &jacobian, NULL, NULL);
-
-  //compute rules
-  decisionRules.compute(jacobian, ghx, ghu);
-}
-void
-ModelSolution::ComputeSteadyState(VectorView &steadyState, const Vector &deepParams)
-{
-  // does nothig for time being.
-}
 
diff --git a/mex/sources/estimation/ModelSolution.hh b/mex/sources/estimation/ModelSolution.hh
index de9773c66..36be5b43a 100644
--- a/mex/sources/estimation/ModelSolution.hh
+++ b/mex/sources/estimation/ModelSolution.hh
@@ -42,7 +42,17 @@ public:
                 const std::vector<size_t> &zeta_back_arg, const std::vector<size_t> &zeta_mixed_arg,
                 const std::vector<size_t> &zeta_static_arg, double qz_criterium);
   virtual ~ModelSolution() {};
-  void compute(VectorView &steadyState, const Vector &deepParams,      Matrix &ghx, Matrix &ghu) throw (DecisionRules::BlanchardKahnException, GeneralizedSchurDecomposition::GSDException);
+  template <class VEC>
+  void compute(VEC &steadyState, const VectorView &deepParams,      Matrix &ghx, Matrix &ghu) throw (DecisionRules::BlanchardKahnException, GeneralizedSchurDecomposition::GSDException)
+  {
+    // compute Steady State
+    ComputeSteadyState(steadyState, deepParams);
+
+    // then get jacobian and
+
+    ComputeModelSolution(steadyState, deepParams, ghx, ghu);
+
+  }
 
 private:
   const size_t n_endo;
@@ -56,8 +66,34 @@ private:
   DynamicModelDLL dynamicDLLp;
   Vector llXsteadyState;
   //Matrix jacobian;
-  void ComputeModelSolution(VectorView &steadyState, const Vector &deepParams,         Matrix &ghx, Matrix &ghu) throw (DecisionRules::BlanchardKahnException, GeneralizedSchurDecomposition::GSDException);
-  void ComputeSteadyState(VectorView &steadyState, const Vector &deepParams);
+  template <class VEC>
+  void ComputeModelSolution(VEC &steadyState, const VectorView &deepParams,         
+			    Matrix &ghx, Matrix &ghu) 
+    throw (DecisionRules::BlanchardKahnException, GeneralizedSchurDecomposition::GSDException)
+  {
+    // set extended Steady State
+
+    for (size_t i = 0; i < zeta_back_mixed.size(); i++)
+      llXsteadyState(i) = steadyState(zeta_back_mixed[i]);
+
+    for (size_t i = 0; i < n_endo; i++)
+      llXsteadyState(zeta_back_mixed.size() + i) = steadyState(i);
+
+    for (size_t i = 0; i < zeta_fwrd_mixed.size(); i++)
+      llXsteadyState(zeta_back_mixed.size() + n_endo + i) = steadyState(zeta_fwrd_mixed[i]);
+
+    //get jacobian
+    dynamicDLLp.eval(llXsteadyState, Mx, deepParams, steadyState, residual, &jacobian, NULL, NULL);
+
+    //compute rules
+    decisionRules.compute(jacobian, ghx, ghu);
+  }
+  template <class VEC>
+  void ComputeSteadyState(VEC &steadyState, const VectorView &deepParams)
+  {
+    // does nothig for time being.
+  }
+
 
 };
 
diff --git a/mex/sources/estimation/RandomWalkMetropolisHastings.cc b/mex/sources/estimation/RandomWalkMetropolisHastings.cc
deleted file mode 100644
index 8e09eb3f4..000000000
--- a/mex/sources/estimation/RandomWalkMetropolisHastings.cc
+++ /dev/null
@@ -1,92 +0,0 @@
-/*
- * Copyright (C) 2010-2011 Dynare Team
- *
- * This file is part of Dynare.
- *
- * Dynare is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * Dynare is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
- */
-
-#include "RandomWalkMetropolisHastings.hh"
-
-#include <iostream>
-#include <fstream>
-
-double
-RandomWalkMetropolisHastings::compute(VectorView &mhLogPostDens, MatrixView &mhParams, Matrix &steadyState,
-                                      Vector &estParams, Vector &deepParams, const MatrixConstView &data, Matrix &Q, Matrix &H,
-                                      const size_t presampleStart, int &info, const size_t startDraw, size_t nMHruns,
-                                      LogPosteriorDensity &lpd, Proposal &pDD, EstimatedParametersDescription &epd)
-{
-  //streambuf *likbuf, *drawbuf *backup;
-  std::ofstream urandfilestr, drawfilestr;
-  urandfilestr.open("urand.csv");
-  drawfilestr.open("paramdraws.csv");
-
-  bool overbound;
-  double newLogpost, logpost, urand;
-  size_t count, accepted = 0;
-  parDraw = estParams;
-
-  logpost = -lpd.compute(steadyState, estParams, deepParams, data, Q, H, presampleStart, info);
-
-  for (size_t run = startDraw - 1; run < nMHruns; ++run)
-    {
-      overbound = false;
-      pDD.draw(parDraw, newParDraw);
-      for (count = 0; count < parDraw.getSize(); ++count)
-        {
-          overbound = (newParDraw(count) < epd.estParams[count].lower_bound || newParDraw(count) > epd.estParams[count].upper_bound);
-          if (overbound)
-            {
-              newLogpost = -INFINITY;
-              break;
-            }
-        }
-      if (!overbound)
-        {
-          try
-            {
-              newLogpost = -lpd.compute(steadyState, newParDraw, deepParams, data, Q, H, presampleStart, info);
-            }
-          catch (const std::exception &e)
-            {
-              throw; // for now handle the system and other errors higher-up
-            }
-          catch (...)
-            {
-              newLogpost = -INFINITY;
-            }
-        }
-      urand = pDD.selectionTestDraw();
-      if ((newLogpost > -INFINITY) && log(urand) < newLogpost-logpost)
-        {
-          parDraw = newParDraw;
-          logpost = newLogpost;
-          accepted++;
-        }
-      mat::get_row(mhParams, run) = parDraw;
-      mhLogPostDens(run) = logpost;
-
-      urandfilestr << urand << "\n"; //","
-      for (size_t c = 0; c < newParDraw.getSize()-1; ++c)
-        drawfilestr << newParDraw(c) << ",";
-      drawfilestr <<  newParDraw(newParDraw.getSize()-1) << "\n";
-    }
-
-  urandfilestr.close();
-  drawfilestr.close();
-
-  return (double) accepted/(nMHruns-startDraw+1);
-}
-
diff --git a/mex/sources/estimation/RandomWalkMetropolisHastings.hh b/mex/sources/estimation/RandomWalkMetropolisHastings.hh
index 5c0578913..38866e76c 100644
--- a/mex/sources/estimation/RandomWalkMetropolisHastings.hh
+++ b/mex/sources/estimation/RandomWalkMetropolisHastings.hh
@@ -26,6 +26,8 @@
 #if !defined(A6BBC5E0_598E_4863_B7FF_E87320056B80__INCLUDED_)
 #define A6BBC5E0_598E_4863_B7FF_E87320056B80__INCLUDED_
 
+#include <iostream>
+#include <fstream>
 #include "LogPosteriorDensity.hh"
 #include "Proposal.hh"
 
@@ -41,11 +43,75 @@ public:
   {
   };
   virtual ~RandomWalkMetropolisHastings() {};
-  virtual double compute(VectorView &mhLogPostDens, MatrixView &mhParams, Matrix &steadyState,
-                         Vector &estParams, Vector &deepParams, const MatrixConstView &data, Matrix &Q, Matrix &H,
-                         const size_t presampleStart, int &info, const size_t startDraw, size_t nMHruns,
-                         LogPosteriorDensity &logPosteriorDensity, Proposal &proposalDrawDistribution,
-                         EstimatedParametersDescription &epd);
+
+  template<class VEC1>
+  double compute(VectorView &mhLogPostDens, MatrixView &mhParams, VEC1 &steadyState,
+		 Vector &estParams, VectorView &deepParams, const MatrixConstView &data, MatrixView &Q, Matrix &H,
+		 const size_t presampleStart, int &info, const size_t startDraw, size_t nMHruns,
+		 LogPosteriorDensity &lpd, Proposal &pDD, EstimatedParametersDescription &epd)
+  {
+    //streambuf *likbuf, *drawbuf *backup;
+    std::ofstream urandfilestr, drawfilestr;
+    urandfilestr.open("urand.csv");
+    drawfilestr.open("paramdraws.csv");
+
+    bool overbound;
+    double newLogpost, logpost, urand;
+    size_t count, accepted = 0;
+    parDraw = estParams;
+
+    logpost = -lpd.compute(steadyState, estParams, deepParams, data, Q, H, presampleStart, info);
+
+    for (size_t run = startDraw - 1; run < nMHruns; ++run)
+      {
+	overbound = false;
+	pDD.draw(parDraw, newParDraw);
+	for (count = 0; count < parDraw.getSize(); ++count)
+	  {
+	    overbound = (newParDraw(count) < epd.estParams[count].lower_bound || newParDraw(count) > epd.estParams[count].upper_bound);
+	    if (overbound)
+	      {
+		newLogpost = -INFINITY;
+		break;
+	      }
+	  }
+	if (!overbound)
+	  {
+	    try
+	      {
+		newLogpost = -lpd.compute(steadyState, newParDraw, deepParams, data, Q, H, presampleStart, info);
+	      }
+	    catch (const std::exception &e)
+	      {
+		throw; // for now handle the system and other errors higher-up
+	      }
+	    catch (...)
+	      {
+		newLogpost = -INFINITY;
+	      }
+	  }
+	urand = pDD.selectionTestDraw();
+	if ((newLogpost > -INFINITY) && log(urand) < newLogpost-logpost)
+	  {
+	    parDraw = newParDraw;
+	    logpost = newLogpost;
+	    accepted++;
+	  }
+	mat::get_row(mhParams, run) = parDraw;
+	mhLogPostDens(run) = logpost;
+
+	urandfilestr << urand << "\n"; //","
+	for (size_t c = 0; c < newParDraw.getSize()-1; ++c)
+	  drawfilestr << newParDraw(c) << ",";
+	drawfilestr <<  newParDraw(newParDraw.getSize()-1) << "\n";
+      }
+
+    urandfilestr.close();
+    drawfilestr.close();
+
+    return (double) accepted/(nMHruns-startDraw+1);
+  };
+
 };
 
 #endif // !defined(A6BBC5E0_598E_4863_B7FF_E87320056B80__INCLUDED_)
diff --git a/mex/sources/estimation/logMHMCMCposterior.cc b/mex/sources/estimation/logMHMCMCposterior.cc
index a77e55c77..c388767af 100644
--- a/mex/sources/estimation/logMHMCMCposterior.cc
+++ b/mex/sources/estimation/logMHMCMCposterior.cc
@@ -114,8 +114,8 @@ fillEstParamsInfo(const mxArray *estim_params_info, EstimatedParameter::pType ty
 
 int
 sampleMHMC(LogPosteriorDensity &lpd, RandomWalkMetropolisHastings &rwmh,
-           Matrix &steadyState, Vector &estParams, Vector &deepParams, const MatrixConstView &data,
-           Matrix &Q, Matrix &H, size_t presampleStart, int &info, const VectorConstView &nruns,
+           VectorView &steadyState, VectorConstView &estParams, VectorView &deepParams, const MatrixConstView &data,
+           MatrixView &Q, Matrix &H, size_t presampleStart, int &info, const VectorConstView &nruns,
            size_t fblock, size_t nBlocks, Proposal pdd, EstimatedParametersDescription &epd,
            const std::string &resultsFileStem, size_t console_mode, size_t load_mh_file)
 {
@@ -539,26 +539,26 @@ sampleMHMC(LogPosteriorDensity &lpd, RandomWalkMetropolisHastings &rwmh,
           mhFName = ssFName.str();
           fidlog = fopen(mhFName.c_str(), "a");
           fprintf(fidlog, "\n");
-          fprintf(fidlog, "%% Mh%dBlck%d ( %s %s )\n", (int) NewFileVw(b-1), b, __DATE__, __TIME__);
+          fprintf(fidlog, "%% Mh%dBlck%lu ( %s %s )\n", (int) NewFileVw(b-1), b, __DATE__, __TIME__);
           fprintf(fidlog, " \n");
-          fprintf(fidlog, "  Number of simulations.: %d \n", currInitSizeArray); // (length(logpo2)) ');
+          fprintf(fidlog, "  Number of simulations.: %lu \n", currInitSizeArray); // (length(logpo2)) ');
           fprintf(fidlog, "  Acceptation rate......: %f \n", jsux);
           fprintf(fidlog, "  Posterior mean........:\n");
           for (size_t i = 0; i < npar; ++i)
             {
               VectorView mhpdColVw = mat::get_col(mhParamDraws, i);
-              fprintf(fidlog, "    params: %d : %f \n", i+1, vec::meanSumMinMax(dsum, dmin, dmax, mhpdColVw));
+              fprintf(fidlog, "    params: %lu : %f \n", i+1, vec::meanSumMinMax(dsum, dmin, dmax, mhpdColVw));
               MinMax(i, iMin) = dmin;
               MinMax(i, iMax) = dmax;
             } // end
           fprintf(fidlog, "    log2po: %f \n", vec::meanSumMinMax(dsum, dmin, dmax, mhLogPostDens));
           fprintf(fidlog, "  Minimum value.........:\n");;
           for (size_t i = 0; i < npar; ++i)
-            fprintf(fidlog, "    params: %d : %f \n", i+1, MinMax(i, iMin));
+            fprintf(fidlog, "    params: %lu : %f \n", i+1, MinMax(i, iMin));
           fprintf(fidlog, "    log2po: %f \n", dmin);
           fprintf(fidlog, "  Maximum value.........:\n");
           for (size_t i = 0; i < npar; ++i)
-            fprintf(fidlog, "    params: %d : %f \n", i+1, MinMax(i, iMax));
+            fprintf(fidlog, "    params: %lu : %f \n", i+1, MinMax(i, iMax));
           fprintf(fidlog, "    log2po: %f \n", dmax);
           fprintf(fidlog, " \n");
           fclose(fidlog);
@@ -624,12 +624,12 @@ sampleMHMC(LogPosteriorDensity &lpd, RandomWalkMetropolisHastings &rwmh,
 }
 
 int
-logMCMCposterior(const VectorConstView &estParams, const MatrixConstView &data,
-                 const size_t fblock, const size_t nBlocks, const VectorConstView &nMHruns, const MatrixConstView &D)
+logMCMCposterior(VectorConstView &estParams, const MatrixConstView &data,
+                 const size_t fblock, const size_t nBlocks, const VectorConstView &nMHruns, const MatrixConstView &D,
+		 VectorView &steadyState, VectorView &deepParams, MatrixView &Q, Matrix &H)
 {
   // Retrieve pointers to global variables
   const mxArray *M_ = mexGetVariablePtr("global", "M_");
-  const mxArray *oo_ = mexGetVariablePtr("global", "oo_");
   const mxArray *options_ = mexGetVariablePtr("global", "options_");
   const mxArray *estim_params_ = mexGetVariablePtr("global", "estim_params_");
 
@@ -643,7 +643,7 @@ logMCMCposterior(const VectorConstView &estParams, const MatrixConstView &data,
 
   size_t n_endo = (size_t) *mxGetPr(mxGetField(M_, 0, "endo_nbr"));
   size_t n_exo = (size_t) *mxGetPr(mxGetField(M_, 0, "exo_nbr"));
-  size_t n_param = (size_t) *mxGetPr(mxGetField(M_, 0, "param_nbr"));
+
   size_t n_estParams = estParams.getSize();
 
   std::vector<size_t> zeta_fwrd, zeta_back, zeta_mixed, zeta_static;
@@ -702,26 +702,8 @@ logMCMCposterior(const VectorConstView &estParams, const MatrixConstView &data,
   LogPosteriorDensity lpd(dynamicDllFile, epd, n_endo, n_exo, zeta_fwrd, zeta_back, zeta_mixed, zeta_static,
                           qz_criterium, varobs, riccati_tol, lyapunov_tol, info);
 
-  // Construct arguments of compute() method
-  Matrix steadyState(n_endo, 1);
-  mat::get_col(steadyState, 0) = VectorConstView(mxGetPr(mxGetField(oo_, 0, "steady_state")), n_endo, 1);
-
-  Vector estParams2(n_estParams);
-  estParams2 = estParams;
-  Vector deepParams(n_param);
-  deepParams = VectorConstView(mxGetPr(mxGetField(M_, 0, "params")), n_param, 1);
-  Matrix Q(n_exo);
-  Q = MatrixConstView(mxGetPr(mxGetField(M_, 0, "Sigma_e")), n_exo, n_exo, n_exo);
-
-  Matrix H(n_varobs);
-  const mxArray *H_mx = mxGetField(M_, 0, "H");
-  if (mxGetM(H_mx) == 1 && mxGetN(H_mx) == 1 && *mxGetPr(H_mx) == 0)
-    H.setAll(0.0);
-  else
-    H = MatrixConstView(mxGetPr(mxGetField(M_, 0, "H")), n_varobs, n_varobs, n_varobs);
-
   // Construct MHMCMC Sampler
-  RandomWalkMetropolisHastings rwmh(estParams2.getSize());
+  RandomWalkMetropolisHastings rwmh(estParams.getSize());
   // Construct GaussianPrior drawDistribution m=0, sd=1
   GaussianPrior drawGaussDist01(0.0, 1.0, -INFINITY, INFINITY, 0.0, 1.0);
   // get Jscale = diag(bayestopt_.jscale);
@@ -731,7 +713,7 @@ logMCMCposterior(const VectorConstView &estParams, const MatrixConstView &data,
   Proposal pdd(vJscale, D);
 
   //sample MHMCMC draws and get get last line run in the last MH block sub-array
-  int lastMHblockArrayLine = sampleMHMC(lpd, rwmh, steadyState, estParams2, deepParams, data, Q, H, presample, info,
+  int lastMHblockArrayLine = sampleMHMC(lpd, rwmh, steadyState, estParams, deepParams, data, Q, H, presample, info,
                                         nMHruns, fblock, nBlocks, pdd, epd, resultsFileStem, console_mode, load_mh_file);
 
   // Cleanups
@@ -746,8 +728,8 @@ void
 mexFunction(int nlhs, mxArray *plhs[],
             int nrhs, const mxArray *prhs[])
 {
-  if (nrhs != 6)
-    DYN_MEX_FUNC_ERR_MSG_TXT("logposterior: exactly six arguments are required.");
+  if (nrhs != 10)
+    DYN_MEX_FUNC_ERR_MSG_TXT("logposterior: exactly 11 arguments are required.");
   if (nlhs != 2)
     DYN_MEX_FUNC_ERR_MSG_TXT("logposterior: exactly two return arguments are required.");
 
@@ -761,18 +743,42 @@ mexFunction(int nlhs, mxArray *plhs[],
   if (!mxIsDouble(prhs[1]))
     DYN_MEX_FUNC_ERR_MSG_TXT("logposterior: Second argument must be a matrix of double-precision numbers");
 
-  MatrixConstView data(mxGetPr(prhs[1]), mxGetM(prhs[1]), mxGetN(prhs[1]), mxGetM(prhs[1]));
-
   size_t fblock = (size_t) mxGetScalar(prhs[2]);
   size_t nBlocks = (size_t) mxGetScalar(prhs[3]);
   VectorConstView nMHruns(mxGetPr(prhs[4]), mxGetM(prhs[4]), 1);
+  MatrixConstView D(mxGetPr(prhs[5]), mxGetM(prhs[5]), mxGetN(prhs[5]), mxGetM(prhs[5]));
+  const mxArray *dataset = prhs[6];
+  const mxArray *options_ = prhs[7];
+  const mxArray *M_ = prhs[8];
+  const mxArray *bayestopt_ = prhs[9];
+  const mxArray *oo_ = prhs[10];
+
   assert(nMHruns.getSize() == nBlocks);
 
-  MatrixConstView D(mxGetPr(prhs[5]), mxGetM(prhs[5]), mxGetN(prhs[5]), mxGetM(prhs[5]));
+  mxArray *dataset_data = mxGetField(dataset,0,"data");
+  MatrixConstView data(mxGetPr(dataset_data), mxGetM(dataset_data), mxGetN(dataset_data), mxGetM(dataset_data));
+
+  int endo_nbr = *(int*)mxGetPr(mxGetField(M_, 0, "endo_nbr"));
+  int exo_nbr = *(int*)mxGetPr(mxGetField(M_, 0, "exo_nbr"));
+  int param_nbr = *(int*)mxGetPr(mxGetField(M_, 0, "param_nbr"));
+  int varobs_nbr = mxGetN(mxGetField(options_, 0, "varobs"));
+
+  VectorView steadyState(mxGetPr(mxGetField(oo_,0,"steady_state")),endo_nbr, 1);
+  VectorView deepParams(mxGetPr(mxGetField(M_, 0, "params")),param_nbr,1);
+
+  MatrixView Q(mxGetPr(mxGetField(M_, 0, "Sigma_e")), exo_nbr, exo_nbr, exo_nbr);
+
+  Matrix H(varobs_nbr,varobs_nbr);
+  const mxArray *H_mx = mxGetField(M_, 0, "H");
+  if (mxGetM(H_mx) == 1 && mxGetN(H_mx) == 1 && *mxGetPr(H_mx) == 0)
+    H.setAll(0.0);
+  else
+    H = MatrixConstView(mxGetPr(H_mx), varobs_nbr, varobs_nbr, varobs_nbr);
+
   //calculate MHMCMC draws and get get last line run in the last MH block sub-array
   try
     {
-      int lastMHblockArrayLine = logMCMCposterior(estParams, data, fblock, nBlocks, nMHruns, D);
+      int lastMHblockArrayLine = logMCMCposterior(estParams, data, fblock, nBlocks, nMHruns, D, steadyState, deepParams, Q, H);
       plhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL);
       *mxGetPr(plhs[1]) = (double) lastMHblockArrayLine;
     }
diff --git a/mex/sources/estimation/logposterior.cc b/mex/sources/estimation/logposterior.cc
index 2f5307647..039a8cc5f 100644
--- a/mex/sources/estimation/logposterior.cc
+++ b/mex/sources/estimation/logposterior.cc
@@ -21,6 +21,7 @@
 #include <vector>
 #include <algorithm>
 #include <functional>
+#include <sstream>
 
 #include "Vector.hh"
 #include "Matrix.hh"
@@ -39,11 +40,10 @@ public:
 };
 
 void
-fillEstParamsInfo(const mxArray *estim_params_info, EstimatedParameter::pType type,
+fillEstParamsInfo(const mxArray *bayestopt_, const mxArray *estim_params_info, EstimatedParameter::pType type,
                   std::vector<EstimatedParameter> &estParamsInfo)
 {
   // execute once only
-  static const mxArray *bayestopt_ = mexGetVariablePtr("global", "bayestopt_");
   static const mxArray *bayestopt_ubp = mxGetField(bayestopt_, 0, "ub"); // upper bound
   static const mxArray *bayestopt_lbp = mxGetField(bayestopt_, 0, "lb"); // lower bound
   static const mxArray *bayestopt_p1p = mxGetField(bayestopt_, 0, "p1"); // prior mean
@@ -100,15 +100,13 @@ fillEstParamsInfo(const mxArray *estim_params_info, EstimatedParameter::pType ty
     }
 }
 
+template <class VEC1, class VEC2>
 double
-logposterior(const VectorConstView &estParams, const MatrixConstView &data)
+logposterior(VEC1 &estParams, const MatrixConstView &data,
+             const mxArray *options_, const mxArray *M_, const mxArray *estim_params_,
+	     const mxArray *bayestopt_, const mxArray *oo_, VEC2 &steadyState, double *trend_coeff,
+	     int &info, VectorView &deepParams, Matrix &H, MatrixView &Q)
 {
-  // Retrieve pointers to global variables
-  const mxArray *M_ = mexGetVariablePtr("global", "M_");
-  const mxArray *oo_ = mexGetVariablePtr("global", "oo_");
-  const mxArray *options_ = mexGetVariablePtr("global", "options_");
-  const mxArray *estim_params_ = mexGetVariablePtr("global", "estim_params_");
-
   // Construct arguments of constructor of LogLikelihoodMain
   char *fName = mxArrayToString(mxGetField(M_, 0, "fname"));
   std::string dynamicDllFile(fName);
@@ -150,6 +148,7 @@ logposterior(const VectorConstView &estParams, const MatrixConstView &data)
     throw LogposteriorMexErrMsgTxtException("options_.varobs_id must be a row vector");
 
   size_t n_varobs = mxGetN(varobs_mx);
+  // substract 1.0 from obsverved variables index
   std::transform(mxGetPr(varobs_mx), mxGetPr(varobs_mx) + n_varobs, back_inserter(varobs),
                  std::bind2nd(std::minus<size_t>(), 1));
 
@@ -160,44 +159,27 @@ logposterior(const VectorConstView &estParams, const MatrixConstView &data)
   estSubsamples.push_back(EstimationSubsample(0, data.getCols() - 1));
 
   std::vector<EstimatedParameter> estParamsInfo;
-  fillEstParamsInfo(mxGetField(estim_params_, 0, "var_exo"), EstimatedParameter::shock_SD,
+  fillEstParamsInfo(bayestopt_, mxGetField(estim_params_, 0, "var_exo"), EstimatedParameter::shock_SD,
                     estParamsInfo);
-  fillEstParamsInfo(mxGetField(estim_params_, 0, "var_endo"), EstimatedParameter::measureErr_SD,
+  fillEstParamsInfo(bayestopt_, mxGetField(estim_params_, 0, "var_endo"), EstimatedParameter::measureErr_SD,
                     estParamsInfo);
-  fillEstParamsInfo(mxGetField(estim_params_, 0, "corrx"), EstimatedParameter::shock_Corr,
+  fillEstParamsInfo(bayestopt_, mxGetField(estim_params_, 0, "corrx"), EstimatedParameter::shock_Corr,
                     estParamsInfo);
-  fillEstParamsInfo(mxGetField(estim_params_, 0, "corrn"), EstimatedParameter::measureErr_Corr,
+  fillEstParamsInfo(bayestopt_, mxGetField(estim_params_, 0, "corrn"), EstimatedParameter::measureErr_Corr,
                     estParamsInfo);
-  fillEstParamsInfo(mxGetField(estim_params_, 0, "param_vals"), EstimatedParameter::deepPar,
+  fillEstParamsInfo(bayestopt_, mxGetField(estim_params_, 0, "param_vals"), EstimatedParameter::deepPar,
                     estParamsInfo);
 
   EstimatedParametersDescription epd(estSubsamples, estParamsInfo);
 
   // Allocate LogPosteriorDensity object
-  int info;
   LogPosteriorDensity lpd(dynamicDllFile, epd, n_endo, n_exo, zeta_fwrd, zeta_back, zeta_mixed, zeta_static,
                           qz_criterium, varobs, riccati_tol, lyapunov_tol, info);
 
   // Construct arguments of compute() method
-  Matrix steadyState(n_endo, 1);
-  mat::get_col(steadyState, 0) = VectorConstView(mxGetPr(mxGetField(oo_, 0, "steady_state")), n_endo, 1);
-
-  Vector estParams2(n_estParams);
-  estParams2 = estParams;
-  Vector deepParams(n_param);
-  deepParams = VectorConstView(mxGetPr(mxGetField(M_, 0, "params")), n_param, 1);
-  Matrix Q(n_exo);
-  Q = MatrixConstView(mxGetPr(mxGetField(M_, 0, "Sigma_e")), n_exo, n_exo, n_exo);
-
-  Matrix H(n_varobs);
-  const mxArray *H_mx = mxGetField(M_, 0, "H");
-  if (mxGetM(H_mx) == 1 && mxGetN(H_mx) == 1 && *mxGetPr(H_mx) == 0)
-    H.setAll(0.0);
-  else
-    H = MatrixConstView(mxGetPr(mxGetField(M_, 0, "H")), n_varobs, n_varobs, n_varobs);
 
   // Compute the posterior
-  double logPD = lpd.compute(steadyState, estParams2, deepParams, data, Q, H, presample, info);
+  double logPD = lpd.compute(steadyState, estParams, deepParams, data, Q, H, presample, info);
 
   // Cleanups
   for (std::vector<EstimatedParameter>::iterator it = estParamsInfo.begin();
@@ -211,31 +193,80 @@ void
 mexFunction(int nlhs, mxArray *plhs[],
             int nrhs, const mxArray *prhs[])
 {
-  if (nrhs != 2 || nlhs != 2)
-    DYN_MEX_FUNC_ERR_MSG_TXT("logposterior: exactly two input arguments and two output arguments are required.");
+  if (nrhs != 7 )
+    DYN_MEX_FUNC_ERR_MSG_TXT("logposterior: exactly 7 input arguments are required.");
 
-  // Check and retrieve the arguments
+  if (nlhs > 9 )
+    DYN_MEX_FUNC_ERR_MSG_TXT("logposterior returns 8 output arguments at the most.");
+
+  // Check and retrieve the RHS arguments
 
   if (!mxIsDouble(prhs[0]) || mxGetN(prhs[0]) != 1)
     DYN_MEX_FUNC_ERR_MSG_TXT("logposterior: First argument must be a column vector of double-precision numbers");
 
   VectorConstView estParams(mxGetPr(prhs[0]), mxGetM(prhs[0]), 1);
 
-  if (!mxIsDouble(prhs[1]))
-    DYN_MEX_FUNC_ERR_MSG_TXT("logposterior: Second argument must be a matrix of double-precision numbers");
+  for (int i = 1; i < 7; ++i)
+    if (!mxIsStruct(prhs[i]))
+      {
+	std::stringstream msg;
+	msg << "logposterior: argument " << i+1 << " must be a Matlab structure";
+	DYN_MEX_FUNC_ERR_MSG_TXT(msg.str().c_str());
+      }
+
+  const mxArray *dataset = prhs[1];
+  const mxArray *options_ = prhs[2];
+  const mxArray *M_ = prhs[3];
+  const mxArray *estim_params_ = prhs[4];
+  const mxArray *bayestopt_ = prhs[5];
+  const mxArray *oo_ = prhs[6];
+
+  mxArray *dataset_data = mxGetField(dataset,0,"data");
+  MatrixConstView data(mxGetPr(dataset_data), mxGetM(dataset_data), mxGetN(dataset_data), mxGetM(dataset_data));
+
+  // Creaete LHS arguments
+
+  size_t endo_nbr = (size_t) *mxGetPr(mxGetField(M_, 0, "endo_nbr"));
+  size_t exo_nbr = (size_t) *mxGetPr(mxGetField(M_, 0, "exo_nbr"));
+  size_t param_nbr = (size_t) *mxGetPr(mxGetField(M_, 0, "param_nbr"));
+  size_t varobs_nbr = mxGetM(mxGetField(options_, 0, "varobs"));
+  plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
+  plhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL);
+  plhs[2] = mxCreateDoubleMatrix(endo_nbr, 1, mxREAL);
+  plhs[3] = mxCreateDoubleMatrix(varobs_nbr, 1, mxREAL);
+  plhs[4] = mxCreateDoubleMatrix(1, 1, mxREAL);
+  plhs[5] = mxCreateDoubleMatrix(param_nbr, 1, mxREAL);
+  plhs[6] = mxCreateDoubleMatrix(varobs_nbr, varobs_nbr, mxREAL);
+  plhs[7] = mxCreateDoubleMatrix(exo_nbr, exo_nbr, mxREAL);
+  double *lik = mxGetPr(plhs[0]);
+  double *exit_flag = mxGetPr(plhs[1]);
+
+  VectorView steadyState(mxGetPr(mxGetField(oo_,0,"steady_state")),endo_nbr, 1);
+  VectorView deepParams(mxGetPr(mxGetField(M_, 0, "params")),param_nbr,1);
+
+  MatrixView Q(mxGetPr(mxGetField(M_, 0, "Sigma_e")), exo_nbr, exo_nbr, exo_nbr);
+
+  Matrix H(varobs_nbr,varobs_nbr);
+  const mxArray *H_mx = mxGetField(M_, 0, "H");
+  if (mxGetM(H_mx) == 1 && mxGetN(H_mx) == 1 && *mxGetPr(H_mx) == 0)
+    H.setAll(0.0);
+  else
+    H = MatrixConstView(mxGetPr(H_mx), varobs_nbr, varobs_nbr, varobs_nbr);
+
+  double *trend_coeff  = mxGetPr(plhs[3]);
+  double *info_mx  = mxGetPr(plhs[4]);
 
-  MatrixConstView data(mxGetPr(prhs[1]), mxGetM(prhs[1]), mxGetN(prhs[1]), mxGetM(prhs[1]));
-  
   // Compute and return the value
   try
     {
-      double lik = logposterior(estParams, data);
-      plhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL);
-      *mxGetPr(plhs[1]) = lik;
+      int info;
+      *lik = logposterior(estParams, data, options_, M_, estim_params_, bayestopt_, oo_,
+				steadyState, trend_coeff, info, deepParams, H, Q);
+      *info_mx = info;
+      *exit_flag = info;
     }
   catch (LogposteriorMexErrMsgTxtException e)
     {
       DYN_MEX_FUNC_ERR_MSG_TXT(e.getErrMsg());
     }
-  plhs[0] = mxCreateDoubleScalar(0);
 }
diff --git a/mex/sources/estimation/utils/dynamic_dll.cc b/mex/sources/estimation/utils/dynamic_dll.cc
index a498396fd..bff1d92a8 100644
--- a/mex/sources/estimation/utils/dynamic_dll.cc
+++ b/mex/sources/estimation/utils/dynamic_dll.cc
@@ -90,10 +90,3 @@ DynamicModelDLL::~DynamicModelDLL()
 #endif
 }
 
-void
-DynamicModelDLL::eval(const Vector &y, const Matrix &x, const Vector &modParams, VectorView &ySteady,
-                      Vector &residual, Matrix *g1, Matrix *g2, Matrix *g3) throw (TSException)
-{
-  Dynamic(y.getData(), x.getData(), 1, modParams.getData(), ySteady.getData(), 0, residual.getData(),
-          g1 == NULL ? NULL : g1->getData(), g2 == NULL ? NULL : g2->getData(), g3 == NULL ? NULL : g3->getData());
-}
diff --git a/mex/sources/estimation/utils/dynamic_dll.hh b/mex/sources/estimation/utils/dynamic_dll.hh
index c3b95a7f0..c1f28e997 100644
--- a/mex/sources/estimation/utils/dynamic_dll.hh
+++ b/mex/sources/estimation/utils/dynamic_dll.hh
@@ -59,4 +59,11 @@ public:
   //! evaluate Dynamic model DLL
   void eval(const Vector &y, const Matrix &x, const Vector &params, VectorView &ySteady,
             Vector &residual, Matrix *g1, Matrix *g2, Matrix *g3) throw (TSException);
+  template<class VEC>
+  void eval(const Vector &y, const Matrix &x, const VectorView &modParams, VEC &ySteady,
+                      Vector &residual, Matrix *g1, Matrix *g2, Matrix *g3) throw (TSException)
+  {
+    Dynamic(y.getData(), x.getData(), 1, modParams.getData(), ySteady.getData(), 0, residual.getData(),
+	    g1 == NULL ? NULL : g1->getData(), g2 == NULL ? NULL : g2->getData(), g3 == NULL ? NULL : g3->getData());
+  };
 };