diff --git a/.gitignore b/.gitignore
index cc77f62c6a5716f6a6c9dbb6e6ee20d6b691b092..3d85d84169d07baf8bdb4689bbfa266099f80b56 100644
--- a/.gitignore
+++ b/.gitignore
@@ -122,6 +122,11 @@ ylwrap
 # Windows
 /windows/dynare-version.nsi
 
+# Estimation DLL tests
+/mex/sources/estimation/tests/test-dr
+/mex/sources/estimation/libmat/tests/test-qr
+/mex/sources/estimation/libmat/tests/test-gsd
+
 # Misc
 !/matlab/swz/c-code/Makefile
 !/mex/sources/kalman/cc/Makefile
diff --git a/Makefile.am b/Makefile.am
index be4207e93aa466801c22b849151c9cb81d97b82b..f5184a861040dcada2f1e5f9ce0d826e15be41c7 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -2,7 +2,7 @@ SUBDIRS = preprocessor doc tests
 
 if HAVE_BLAS
 if HAVE_LAPACK
-SUBDIRS += dynare++
+SUBDIRS += dynare++ mex/sources/estimation/tests mex/sources/estimation/libmat/tests
 endif
 endif
 
diff --git a/configure.ac b/configure.ac
index e2317beacdcd3fcee43c4c28b40334f1149726a3..f98dd20aae22aef9b83dcf89b6dc5c1b77e5bbfe 100644
--- a/configure.ac
+++ b/configure.ac
@@ -171,6 +171,8 @@ AC_CONFIG_FILES([Makefile
                  dynare++/integ/testing/Makefile
                  dynare++/kord/Makefile
                  dynare++/src/Makefile
+                 mex/sources/estimation/tests/Makefile
+                 mex/sources/estimation/libmat/tests/Makefile
 ])
 
 AC_ARG_ENABLE([matlab], AS_HELP_STRING([--disable-matlab], [disable compilation of MEX files for MATLAB]), [], [enable_matlab=yes])
diff --git a/mex/sources/dynblas.h b/mex/sources/dynblas.h
index 15820e421240dec43a38e55cc6124e528e862746..aea257b6c7ace09232baa149024a16fcfd54d06f 100644
--- a/mex/sources/dynblas.h
+++ b/mex/sources/dynblas.h
@@ -8,7 +8,7 @@
  * and MATLAB_VERSION (for version 7.4, define it to 0x0704).
  *
  *
- * Copyright (C) 2009 Dynare Team
+ * Copyright (C) 2009-2010 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -62,6 +62,12 @@ extern "C" {
              CONST_BLDOU b, CONST_BLINT ldb, CONST_BLDOU beta,
              BLDOU c, CONST_BLINT ldc);
 
+#define dsymm FORTRAN_WRAPPER(dsymm)
+  void dsymm(BLCHAR side, BLCHAR uplo, CONST_BLINT m, CONST_BLINT n,
+             CONST_BLDOU alpha, CONST_BLDOU a, CONST_BLINT lda,
+             CONST_BLDOU b, CONST_BLINT ldb, CONST_BLDOU beta,
+             BLDOU c, CONST_BLINT ldc);
+
 #define dgemv FORTRAN_WRAPPER(dgemv)
   void dgemv(BLCHAR trans, CONST_BLINT m, CONST_BLINT n, CONST_BLDOU alpha,
              CONST_BLDOU a, CONST_BLINT lda, CONST_BLDOU x, CONST_BLINT incx,
@@ -99,6 +105,10 @@ extern "C" {
   double ddot(CONST_BLINT n, CONST_BLDOU x, CONST_BLINT incx, CONST_BLDOU y,
               CONST_BLINT incy);
 
+#define dsyr FORTRAN_WRAPPER(dsyr)
+  void dsyr(BLCHAR uplo, CONST_BLINT n, CONST_BLDOU alpha, CONST_BLDOU x,
+            CONST_BLINT incx, BLDOU a, CONST_BLINT lda);
+
 #ifdef __cplusplus
 } /* extern "C" */
 #endif
diff --git a/mex/sources/dynlapack.h b/mex/sources/dynlapack.h
index 7bed4b31470538ec9d6910fe66b2d4375e12c573..e918ff37e8e734005cab007bd57b8436bfd8cf9a 100644
--- a/mex/sources/dynlapack.h
+++ b/mex/sources/dynlapack.h
@@ -8,7 +8,7 @@
  * and MATLAB_VERSION (for version 7.4, define it to 0x0704).
  *
  *
- * Copyright (C) 2009 Dynare Team
+ * Copyright (C) 2009-2010 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -110,6 +110,15 @@ extern "C" {
               LAINT isuppz, LADOU work, CONST_LAINT lwork, LAINT iwork, CONST_LAINT liwork,
               LAINT info);
 
+#define dgeqrf FORTRAN_WRAPPER(dgeqrf)
+  void dgeqrf(CONST_LAINT m, CONST_LAINT n, LADOU a, CONST_LAINT lda,
+              LADOU tau, LADOU work, CONST_LAINT lwork, LAINT info);
+
+#define dormqr FORTRAN_WRAPPER(dormqr)
+  void dormqr(LACHAR side, LACHAR trans, CONST_LAINT m, CONST_LAINT n, CONST_LAINT k,
+              CONST_LADOU a, CONST_LAINT lda, CONST_LADOU tau, LADOU c, CONST_LAINT ldc,
+              LADOU work, CONST_LAINT lwork, LAINT info);
+
 #ifdef __cplusplus
 } /* extern "C" */
 #endif
diff --git a/mex/sources/estimation/DecisionRules.cc b/mex/sources/estimation/DecisionRules.cc
new file mode 100644
index 0000000000000000000000000000000000000000..9c918c01ceb7fadbd69bbec632839aa196f79607
--- /dev/null
+++ b/mex/sources/estimation/DecisionRules.cc
@@ -0,0 +1,195 @@
+/*
+ * Copyright (C) 2010 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 <cassert>
+
+#include <algorithm>
+
+#include "DecisionRules.hh"
+
+DecisionRules::DecisionRules(size_t n_arg, size_t p_arg,
+                             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,
+                             double qz_criterium) :
+  n(n_arg), p(p_arg), zeta_fwrd(zeta_fwrd_arg), zeta_back(zeta_back_arg),
+  zeta_mixed(zeta_mixed_arg), zeta_static(zeta_static_arg),
+  n_fwrd(zeta_fwrd.size()), n_back(zeta_back.size()),
+  n_mixed(zeta_mixed.size()), n_static(zeta_static.size()),
+  n_back_mixed(n_back+n_mixed), n_fwrd_mixed(n_fwrd+n_mixed),
+  n_dynamic(n-n_static),
+  S(n, n_static),
+  A(n_back_mixed + n + n_fwrd_mixed),
+  D(n_fwrd + n_back + 2*n_mixed),
+  E(n_fwrd + n_back + 2*n_mixed),
+  Z(n_fwrd + n_back + 2*n_mixed),
+  QR(n, n_static),
+  GSD(n_fwrd + n_back + 2*n_mixed, qz_criterium),
+  LU1(n_fwrd_mixed),
+  LU2(n_back_mixed),
+  LU3(n_static),
+  g_y_back(n_back_mixed),
+  g_y_back_tmp(n_back_mixed),
+  g_y_static(n_static, n_back_mixed),
+  A0s(n_static),
+  A0d(n_static, n_dynamic),
+  g_y_dynamic(n_dynamic, n_back_mixed),
+  g_y_static_tmp(n_static, n_back_mixed)
+{
+  assert(n == n_back + n_fwrd + n_mixed + n_static);
+
+  set_union(zeta_fwrd.begin(), zeta_fwrd.end(),
+            zeta_mixed.begin(), zeta_mixed.end(),
+            back_inserter(zeta_fwrd_mixed));
+  set_union(zeta_back.begin(), zeta_back.end(),
+            zeta_mixed.begin(), zeta_mixed.end(),
+            back_inserter(zeta_back_mixed));
+  set_union(zeta_back_mixed.begin(), zeta_back_mixed.end(),
+            zeta_fwrd.begin(), zeta_fwrd.end(),
+            back_inserter(zeta_dynamic));
+
+  // Compute beta_back and pi_back
+  for(size_t i = 0; i < n_back_mixed; i++)
+    if (find(zeta_mixed.begin(), zeta_mixed.end(), zeta_back_mixed[i])
+        == zeta_mixed.end())
+      beta_back.push_back(i);
+    else
+      pi_back.push_back(i);
+
+  // Compute beta_fwrd and pi_fwrd
+  for(size_t i = 0; i < n_fwrd_mixed; i++)
+    if (find(zeta_mixed.begin(), zeta_mixed.end(), zeta_fwrd_mixed[i])
+        == zeta_mixed.end())
+      beta_fwrd.push_back(i);
+    else
+      pi_fwrd.push_back(i);
+}
+
+void
+DecisionRules::compute(const Matrix &jacobian, Matrix &g_y, Matrix &g_u) throw (BlanchardKahnException, GeneralizedSchurDecomposition::GSDException)
+{
+  assert(jacobian.getRows() == n
+         && jacobian.getCols() == (n_back_mixed + n + n_fwrd_mixed + p));
+  assert(g_y.getRows() == n && g_y.getCols() == n);
+  assert(g_u.getRows() == n && g_u.getCols() == p);
+
+  // Construct S, perform QR decomposition and get A = Q*jacobian
+  for(size_t i = 0; i < n_static; i++)
+    mat::col_copy(jacobian, n_back_mixed+zeta_static[i], S, i);
+
+  A = MatrixConstView(jacobian, 0, 0, n, n_back_mixed + n + n_fwrd_mixed);
+  QR.computeAndMultByQ(S, "L", "N", A);
+
+  // Construct matrix D
+  for (size_t j = 0; j < n_back_mixed; j++)
+    mat::col_copy(A, n_back_mixed + zeta_back_mixed[j], n_static, n - n_static,
+                  D, j, 0);
+  MatrixView(D, 0, n_back_mixed, n - n_static, n_fwrd_mixed) = MatrixView(A, n_static, n_back_mixed + n, n - n_static, n_fwrd_mixed);
+  MatrixView(D, n - n_static, 0, n_mixed, n_fwrd + n_back + 2*n_mixed).setAll(0.0);
+  for (size_t i = 0; i < n_mixed; i++)
+    D(n - n_static + i, beta_back[i]) = 1.0;
+
+  // Construct matrix E
+  MatrixView(E, 0, 0, n - n_static, n_back_mixed) = MatrixView(A, n_static, 0, n - n_static, n_back_mixed);
+  for (size_t j = 0; j < n_fwrd_mixed; j++)
+    mat::col_copy(A, n_back_mixed + zeta_fwrd_mixed[j], n_static, n - n_static,
+                  E, n_back_mixed + j, 0);
+  mat::negate(MatrixView(E, 0, 0, n - n_static, n_fwrd + n_back + 2*n_mixed));
+  MatrixView(E, n - n_static, 0, n_mixed, n_fwrd + n_back + 2*n_mixed).setAll(0.0);
+  for (size_t i = 0; i < n_mixed; i++)
+    E(n - n_static + i, n_back_mixed + beta_fwrd[i]) = 1.0;
+
+  // Perform the generalized Schur
+  size_t sdim;
+  GSD.compute(D, E, Z, sdim);
+
+  if (n_back_mixed != sdim)
+    throw BlanchardKahnException(true, n_fwrd_mixed, n_fwrd + n_back + 2*n_mixed - sdim);
+
+  // Compute DR for forward variables w.r. to endogenous
+  MatrixView Z21(Z, n_back_mixed, 0, n_fwrd_mixed, n_back_mixed),
+    Z22(Z, n_back_mixed, n_back_mixed, n_fwrd_mixed, n_fwrd_mixed);
+
+  try
+    {
+      LU1.invMult("N", Z22, Z21);
+    }
+  catch (LUSolver::LUException &e)
+    {
+      throw BlanchardKahnException(false, n_fwrd_mixed, n_fwrd + n_back + 2*n_mixed - sdim);
+    }
+  mat::negate(Z21);
+  const MatrixView &g_y_fwrd = Z21;
+
+  for (size_t i = 0; i < n_fwrd_mixed; i++)
+    mat::row_copy(g_y_fwrd, i, g_y, zeta_fwrd_mixed[i]);
+
+  // Compute DR for backward variables w.r. to endogenous
+  MatrixView Z11(Z, 0, 0, n_back_mixed, n_back_mixed),
+    T11(D, 0, 0, n_back_mixed, n_back_mixed),
+    S11(E, 0, 0, n_back_mixed, n_back_mixed);
+  mat::set_identity(g_y_back);
+  mat::transpose(Z11);
+  g_y_back_tmp = Z11;
+  LU2.invMult("N", g_y_back_tmp, g_y_back);
+  g_y_back_tmp = g_y_back;
+  blas::gemm("N", "N", 1.0, S11, g_y_back_tmp, 0.0, g_y_back);
+  LU2.invMult("N", T11, g_y_back);
+  g_y_back_tmp = g_y_back;
+  blas::gemm("N", "N", 1.0, Z11, g_y_back_tmp, 0.0, g_y_back);
+
+  // TODO: avoid to copy mixed variables again, rather test it...
+  for (size_t i = 0; i < n_back_mixed; i++)
+    mat::row_copy(g_y_back, i, g_y, zeta_back_mixed[i]);
+
+  // Compute DR for static variables w.r. to endogenous
+  g_y_static = MatrixView(A, 0, 0, n_static, n_back_mixed);
+  for (size_t i = 0; i < n_dynamic; i++)
+    {
+      mat::row_copy(g_y, zeta_dynamic[i], g_y_dynamic, i);
+      mat::col_copy(A, n_back_mixed + zeta_dynamic[i], A0d, i);
+    }
+  blas::gemm("N", "N", 1.0, A0d, g_y_dynamic, 1.0, g_y_static);
+  blas::gemm("N", "N", 1.0, g_y_fwrd, g_y_back, 0.0, g_y_static_tmp);
+  blas::gemm("N", "N", 1.0, MatrixView(A, 0, n_back_mixed + n, n_static, n_fwrd_mixed),
+             g_y_static_tmp, 1.0, g_y_static);
+  for (size_t i = 0; i < n_static; i++)
+    mat::col_copy(A, n_back_mixed + zeta_static[i], A0s, i);
+  LU3.invMult("N", A0s, g_y_static);
+  mat::negate(g_y_static);
+
+  for (size_t i = 0; i < n_static; i++)
+    mat::row_copy(g_y_static, i, g_y, zeta_static[i]);
+
+  // Compute DR for all endogenous w.r. to shocks
+
+}
+
+std::ostream &
+operator<<(std::ostream &out, const DecisionRules::BlanchardKahnException &e)
+{
+  if (e.order)
+    out << "The Blanchard-Kahn order condition is not satisfied: you have " << e.n_fwrd_vars << " forward variables for " << e.n_explosive_eigenvals << " explosive eigenvalues";
+  else
+    out << "The Blanchard Kahn rank condition is not satisfied";
+  out << std::endl;
+  return out;
+}
+
diff --git a/mex/sources/estimation/DecisionRules.hh b/mex/sources/estimation/DecisionRules.hh
new file mode 100644
index 0000000000000000000000000000000000000000..e81da8bfc82207ce0d6884e0d1d61f0315a1f4cb
--- /dev/null
+++ b/mex/sources/estimation/DecisionRules.hh
@@ -0,0 +1,64 @@
+/*
+ * Copyright (C) 2010 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 <cstdlib>
+#include <vector>
+
+#include "Vector.hh"
+#include "Matrix.hh"
+#include "QRDecomposition.hh"
+#include "GeneralizedSchurDecomposition.hh"
+#include "LUSolver.hh"
+
+class DecisionRules
+{
+private:
+  const size_t n, p;
+  const std::vector<size_t> zeta_fwrd, zeta_back, zeta_mixed, zeta_static;
+  const size_t n_fwrd, n_back, n_mixed, n_static, n_back_mixed, n_fwrd_mixed, n_dynamic;
+  std::vector<size_t> zeta_fwrd_mixed, zeta_back_mixed, zeta_dynamic,
+    beta_back, beta_fwrd, pi_back, pi_fwrd;
+  Matrix S, A, D, E, Z;
+  QRDecomposition QR;
+  GeneralizedSchurDecomposition GSD;
+  LUSolver LU1, LU2, LU3;
+  Matrix g_y_back, g_y_back_tmp;
+  Matrix g_y_static, A0s, A0d, g_y_dynamic, g_y_static_tmp;
+public:
+  class BlanchardKahnException
+  {
+  public:
+    //! True if the model fails the order condition. False if it fails the rank condition.
+    const bool order;
+    const int n_fwrd_vars, n_explosive_eigenvals;
+    BlanchardKahnException(bool order_arg, int n_fwrd_vars_arg, int n_explosive_eigenvals_arg) : order(order_arg), n_fwrd_vars(n_fwrd_vars_arg), n_explosive_eigenvals(n_explosive_eigenvals_arg) {};
+  };
+  /*!
+    The zetas are supposed to follow C convention (first vector index is zero).
+  */
+  DecisionRules(size_t n_arg, size_t p_arg, 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, double qz_criterium);
+  /*!
+    \param jacobian First columns are backetermined vars at t-1 (in the order of zeta_back_mixed), then all vars at t (in the orig order), then forward vars at t+1 (in the order of zeta_fwrd_mixed), then exogenous vars.
+  */
+  void compute(const Matrix &jacobian, Matrix &g_y, Matrix &g_u) throw (BlanchardKahnException, GeneralizedSchurDecomposition::GSDException);
+};
+
+std::ostream &operator<<(std::ostream &out, const DecisionRules::BlanchardKahnException &e);
diff --git a/mex/sources/estimation/libmat/BlasBindings.hh b/mex/sources/estimation/libmat/BlasBindings.hh
new file mode 100644
index 0000000000000000000000000000000000000000..4196b2f88d227fabe63b105ed28c865f58300ef3
--- /dev/null
+++ b/mex/sources/estimation/libmat/BlasBindings.hh
@@ -0,0 +1,80 @@
+/*
+ * Copyright (C) 2010 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/>.
+ */
+
+#ifndef _BLAS_BINDINGS_HH
+#define _BLAS_BINDINGS_HH
+
+#include <dynblas.h>
+
+#include "Vector.hh"
+#include "Matrix.hh"
+
+namespace blas
+{
+  /* Level 2 */
+  
+  //! Symmetric rank 1 operation: A = alpha*X*X' + A
+  template<class Mat, class Vec>
+  inline void
+  syr(const char *uplo, double alpha, Vec X, Mat A)
+  {
+    assert(X.getSize() == A.getCols() && A.getCols() == A.getRows());
+    blas_int n = X.getSize();
+    blas_int incx = X.getStride();
+    blas_int lda = A.getLd();
+    dsyr(uplo, &n, &alpha, X.getData(), &incx, A.getData(), &lda);
+  }
+
+  /* Level 3 */
+
+  //! General matrix multiplication
+  template<class Mat1, class Mat2, class Mat3>
+  inline void
+  gemm(const char *transa, const char *transb,
+       double alpha, const Mat1 &A, const Mat2 &B,
+       double beta, Mat3 &C)
+  {
+    assert(A.getRows() == C.getRows());
+    assert(A.getCols() == B.getRows());
+    assert(B.getCols() == C.getCols());
+    blas_int m = A.getRows(), n = B.getCols(), k = A.getCols();
+    blas_int lda = A.getLd(), ldb = B.getLd(), ldc = C.getLd();
+    dgemm(transa, transb, &m, &n, &k, &alpha, A.getData(), &lda,
+          B.getData(), &ldb, &beta, C.getData(), &ldc);
+  }
+
+  //! Symmetric matrix multiplication
+  template<class Mat1, class Mat2, class Mat3>
+  inline void
+  symm(const char *side, const char *uplo,
+       double alpha, const Mat1 &A, const Mat2 &B,
+       double beta, Mat3 &C)
+  {
+    assert(A.getRows() == A.getCols());
+    assert(A.getRows() == C.getRows());
+    assert(A.getCols() == B.getRows());
+    assert(B.getCols() == C.getCols());
+    blas_int m = A.getRows(), n = B.getCols();
+    blas_int lda = A.getLd(), ldb = B.getLd(), ldc = C.getLd();
+    dsymm(side, uplo, &m, &n, &alpha, A.getData(), &lda,
+          B.getData(), &ldb, &beta, C.getData(), &ldc);
+  }
+} // End of namespace
+
+#endif
diff --git a/mex/sources/estimation/libmat/GeneralizedSchurDecomposition.cc b/mex/sources/estimation/libmat/GeneralizedSchurDecomposition.cc
new file mode 100644
index 0000000000000000000000000000000000000000..cbd75c534ca5168a2b61c919e16ca3d5bc754bff
--- /dev/null
+++ b/mex/sources/estimation/libmat/GeneralizedSchurDecomposition.cc
@@ -0,0 +1,73 @@
+/*
+ * Copyright (C) 2010 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 "GeneralizedSchurDecomposition.hh"
+
+#include <cassert>
+#include <cstdlib>
+
+double GeneralizedSchurDecomposition::criterium_static;
+
+GeneralizedSchurDecomposition::GeneralizedSchurDecomposition(size_t n_arg, double criterium_arg)
+  : n(n_arg), criterium(criterium_arg)
+{
+  alphar = new double[n];
+  alphai = new double[n];
+  beta = new double[n];
+
+  lwork = 16*n+16; // Same heuristic choice than mjdgges
+  work = new double[(int) lwork];
+  vsl = new double[n*n];
+
+  bwork = new lapack_int[n];
+}
+
+GeneralizedSchurDecomposition::~GeneralizedSchurDecomposition()
+{
+  delete[] alphar;
+  delete[] alphai;
+  delete[] beta;
+  delete[] work;
+  delete[] vsl;
+  delete[] bwork;
+}
+
+lapack_int
+GeneralizedSchurDecomposition::selctg(const double *alphar, const double *alphai, const double *beta)
+{
+  return((*alphar * *alphar + *alphai * *alphai) < criterium_static * *beta * *beta);
+}
+
+std::ostream &
+operator<<(std::ostream &out, const GeneralizedSchurDecomposition::GSDException &e)
+{
+  out << "DGGES return code " << e.info << ": ";
+  if (e.info < 0)
+    out << "argument " << -e.info << " has an illegal value";
+  else if (e.info <= e.n)
+    out << "the QZ iteration failed";
+  else if (e.info == e.n + 1)
+    out << "other than QZ iteration failed in DHGEQZ";
+  else if (e.info == e.n + 2)
+    out << "after reordering, roundoff changed values of some complex eigenvalues so that leading eigenvalues in the Generalized Schur form no longer satisfy SELCTG=TRUE. This could also be caused due to scaling";
+  else if (e.info == e.n + 3)
+    out << "reordering failed in DTGSEN";
+  out << std::endl;
+  return out;
+}
diff --git a/mex/sources/estimation/libmat/GeneralizedSchurDecomposition.hh b/mex/sources/estimation/libmat/GeneralizedSchurDecomposition.hh
new file mode 100644
index 0000000000000000000000000000000000000000..bf366597a8b3ba526761a9747fc83f03e830f078
--- /dev/null
+++ b/mex/sources/estimation/libmat/GeneralizedSchurDecomposition.hh
@@ -0,0 +1,110 @@
+/*
+ * Copyright (C) 2010 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 <dynlapack.h>
+
+#include "Vector.hh"
+#include "Matrix.hh"
+
+class GeneralizedSchurDecomposition
+{
+private:
+  const size_t n;
+  const double criterium;
+  lapack_int lwork;
+  double *alphar, *alphai, *beta, *vsl, *work;
+  lapack_int *bwork;
+  static double criterium_static;
+  static lapack_int selctg(const double *alphar, const double *alphai, const double *beta);
+public:
+  class GSDException
+  {
+  public:
+    const lapack_int info, n;
+    GSDException(lapack_int info_arg, lapack_int n_arg) : info(info_arg), n(n_arg) {};
+  };
+  //! \todo Replace heuristic choice for workspace size by a query to determine the optimal size
+  GeneralizedSchurDecomposition(size_t n_arg, double criterium_arg);
+  ~GeneralizedSchurDecomposition();
+  //! \todo Add a lock around the modification of criterium_static for making it thread-safe
+  template<class Mat1, class Mat2, class Mat3>
+  void compute(Mat1 &S, Mat2 &T, Mat3 &Z, size_t &sdim) throw (GSDException);
+  template<class Mat1, class Mat2, class Mat3, class Mat4, class Mat5>
+  /*!
+    \param[out] sdim Number of non-explosive generalized eigenvalues
+  */
+  void compute(const Mat1 &D, const Mat2 &E, Mat3 &S, Mat4 &T, Mat5 &Z, size_t &sdim) throw (GSDException);
+  template<class Vec1, class Vec2>
+  void getGeneralizedEigenvalues(Vec1 &eig_real, Vec2 &eig_cmplx);
+};
+
+std::ostream &operator<<(std::ostream &out, const GeneralizedSchurDecomposition::GSDException &e);
+
+template<class Mat1, class Mat2, class Mat3>
+void
+GeneralizedSchurDecomposition::compute(Mat1 &S, Mat2 &T, Mat3 &Z, size_t &sdim) throw (GSDException)
+{
+  assert(S.getRows() == n && S.getCols() == n
+         && T.getRows() == n && T.getCols() == n
+         && Z.getRows() == n && Z.getCols() == n);
+
+  lapack_int n2 = n;
+  lapack_int info, sdim2;
+  lapack_int lds = S.getLd(), ldt = T.getLd(), ldz = Z.getLd();
+
+  criterium_static = criterium;
+  // Here we are forced to give space for left Schur vectors, even if we don't use them, because of a bug in dgges()
+  dgges("N", "V", "S", &selctg, &n2, S.getData(), &lds, T.getData(), &ldt,
+        &sdim2, alphar, alphai, beta, vsl, &n2, Z.getData(), &ldz,
+        work, &lwork, bwork, &info);
+
+  if (info != 0)
+    throw GSDException(info, n2);
+
+  sdim = sdim2;
+}
+
+template<class Mat1, class Mat2, class Mat3, class Mat4, class Mat5>
+void
+GeneralizedSchurDecomposition::compute(const Mat1 &D, const Mat2 &E,
+                                       Mat3 &S, Mat4 &T, Mat5 &Z, size_t &sdim) throw (GSDException)
+{
+  assert(D.getRows() == n && D.getCols() == n
+         && E.getRows() == n && E.getCols() == n);
+  S = D;
+  T = E;
+  compute(S, T, Z, sdim);
+}
+
+template<class Vec1, class Vec2>
+void
+GeneralizedSchurDecomposition::getGeneralizedEigenvalues(Vec1 &eig_real, Vec2 &eig_cmplx)
+{
+  assert(eig_real.getSize() == n && eig_cmplx.getSize() == n);
+
+  double *par = alphar, *pai = alphai, *pb = beta,
+    *per = eig_real.getData(), *pei = eig_cmplx.getData();
+  while (par < alphar + n)
+    {
+      *per = *par++ / *pb;
+      *pei = *pai++ / *pb++;
+      per += eig_real.getStride();
+      pei += eig_cmplx.getStride();
+    }
+}
diff --git a/mex/sources/estimation/libmat/LUSolver.cc b/mex/sources/estimation/libmat/LUSolver.cc
new file mode 100644
index 0000000000000000000000000000000000000000..30033f1ebf9bd74dd005d7dd2356b678e9b04a1e
--- /dev/null
+++ b/mex/sources/estimation/libmat/LUSolver.cc
@@ -0,0 +1,30 @@
+/*
+ * Copyright (C) 2010 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 "LUSolver.hh"
+
+LUSolver::LUSolver(size_t dim_arg) : dim(dim_arg)
+{
+  ipiv = new lapack_int[dim];
+}
+
+LUSolver::~LUSolver()
+{
+  delete[] ipiv;
+}
diff --git a/mex/sources/estimation/libmat/LUSolver.hh b/mex/sources/estimation/libmat/LUSolver.hh
new file mode 100644
index 0000000000000000000000000000000000000000..d6d24db56d9ab19a9660148050e0c309c73b7124
--- /dev/null
+++ b/mex/sources/estimation/libmat/LUSolver.hh
@@ -0,0 +1,63 @@
+/*
+ * Copyright (C) 2010 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 <cstdlib>
+#include <cassert>
+
+#include <dynlapack.h>
+
+class LUSolver
+{
+private:
+  const size_t dim;
+  lapack_int *ipiv;
+public:
+  class LUException
+  {
+  public:
+    const lapack_int info;
+    LUException(lapack_int info_arg) : info(info_arg) {};
+  };
+  LUSolver(size_t dim_arg);
+  ~LUSolver();
+  /*!
+    Computes A^(-1)*B (possibly transposing A).
+    The output is stored in B.
+    The input matrix A is modified.
+  */
+  template<class Mat1, class Mat2>
+  void invMult(const char *trans, Mat1 &A, Mat2 &B) throw (LUException);
+};
+
+template<class Mat1, class Mat2>
+void
+LUSolver::invMult(const char *trans, Mat1 &A, Mat2 &B) throw (LUException)
+{
+  assert(A.getRows() == dim && A.getCols() == dim);
+  assert(B.getRows() == dim);
+  lapack_int n = dim, lda = A.getLd(), info;
+  dgetrf(&n, &n, A.getData(), &lda, ipiv, &info);
+
+  if (info != 0)
+    throw LUException(info);
+
+  lapack_int nrhs = B.getCols(), ldb = B.getLd();
+  dgetrs(trans, &n, &nrhs, A.getData(), &lda, ipiv, B.getData(), &ldb, &info);
+  assert(info == 0);
+}
diff --git a/mex/sources/estimation/libmat/Matrix.cc b/mex/sources/estimation/libmat/Matrix.cc
new file mode 100644
index 0000000000000000000000000000000000000000..bc43da795d2b0c9952c73f8ef830f6d1a15065e5
--- /dev/null
+++ b/mex/sources/estimation/libmat/Matrix.cc
@@ -0,0 +1,112 @@
+/*
+ * Copyright (C) 2010 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 "Matrix.hh"
+
+#include <cstring>
+#include <cassert>
+
+Matrix::Matrix(size_t rows_arg, size_t cols_arg) : rows(rows_arg), cols(cols_arg)
+{
+  data = new double[rows*cols];
+}
+
+Matrix::Matrix(size_t size_arg) : rows(size_arg), cols(size_arg)
+{
+  data = new double[rows*cols];
+}
+
+Matrix::Matrix(const Matrix &arg) : rows(arg.rows), cols(arg.cols)
+{
+  data = new double[rows*cols];
+  memcpy(data, arg.data, rows*cols*sizeof(double));
+}
+
+Matrix::~Matrix()
+{
+  delete[] data;
+}
+
+Matrix &
+Matrix::operator= (const Matrix &arg)
+{
+  assert(rows == arg.rows && cols == arg.cols);
+  memcpy(data, arg.data, rows*cols*sizeof(double));
+  return *this;
+}
+
+std::ostream &
+operator<<(std::ostream &out, const Matrix &M)
+{
+  mat::print(out, M);
+  return out;
+}
+
+MatrixView::MatrixView(double *data_arg, size_t rows_arg, size_t cols_arg, size_t ld_arg)
+  : data(data_arg), rows(rows_arg), cols(cols_arg), ld(ld_arg)
+{
+}
+
+MatrixView::MatrixView(Matrix &arg, size_t row_offset, size_t col_offset,
+                       size_t rows_arg, size_t cols_arg) :
+  data(arg.getData() + row_offset + col_offset*arg.getLd()), rows(rows_arg), cols(cols_arg), ld(arg.getLd())
+{
+  assert(row_offset < arg.getRows()
+         && row_offset + rows_arg <= arg.getRows()
+         && col_offset < arg.getCols()
+         && col_offset + cols_arg <= arg.getCols());
+}
+
+std::ostream &
+operator<<(std::ostream &out, const MatrixView &M)
+{
+  mat::print(out, M);
+  return out;
+}
+
+MatrixView &
+MatrixView::operator= (const MatrixView &arg)
+{
+  assert(rows == arg.getRows() && cols == arg.getCols());
+  for (size_t j = 0; j < cols; j++)
+    memcpy(data + j*ld, arg.getData() + j*arg.getLd(), rows*sizeof(double));
+  return *this;
+}
+
+MatrixConstView::MatrixConstView(const double *data_arg, size_t rows_arg, size_t cols_arg, size_t ld_arg)
+  : data(data_arg), rows(rows_arg), cols(cols_arg), ld(ld_arg)
+{
+}
+
+MatrixConstView::MatrixConstView(const Matrix &arg, size_t row_offset, size_t col_offset,
+                                 size_t rows_arg, size_t cols_arg) :
+  data(arg.getData() + row_offset + col_offset*arg.getLd()), rows(rows_arg), cols(cols_arg), ld(arg.getLd())
+{
+  assert(row_offset < arg.getRows()
+         && row_offset + rows_arg <= arg.getRows()
+         && col_offset < arg.getCols()
+         && col_offset + cols_arg <= arg.getCols());
+}
+
+std::ostream &
+operator<<(std::ostream &out, const MatrixConstView &M)
+{
+  mat::print(out, M);
+  return out;
+}
diff --git a/mex/sources/estimation/libmat/Matrix.hh b/mex/sources/estimation/libmat/Matrix.hh
new file mode 100644
index 0000000000000000000000000000000000000000..b1feab6c914ac701badb6e4004017ea9b8954944
--- /dev/null
+++ b/mex/sources/estimation/libmat/Matrix.hh
@@ -0,0 +1,348 @@
+/*
+ * Copyright (C) 2010 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/>.
+ */
+
+#ifndef _MATRIX_HH
+#define _MATRIX_HH
+
+#include <algorithm>
+#include <iostream>
+#include <iomanip>
+
+#include <cstdlib>
+#include <cassert>
+#include <cstring>
+
+#include "Vector.hh"
+
+/*
+  This header defines three matrix classes, which implement a "matrix concept"
+  (much like the concepts of the Standard Template Library or of Boost
+  Library). The first class is a matrix owning the data space for its
+  elements, and the other two are matrix "views" of another matrix, i.e. a
+  contiguous submatrix. This design philosophy is close to the design of the
+  GNU Scientific Library, but here using the syntactic power of C++ class and
+  templates, while achieving very high efficiency.
+
+  These classes can be used with various templated functions, including
+  wrappers around the BLAS primitives.
+
+  The expressions required to be valid for a class M implementing the "matrix concept" are:
+  - M.getRows(): return number of rows
+  - M.getCols(): return number of columns
+  - M.getLd(): return the leading dimension (here the offset between two columns in the data space, since we are in column-major order)
+  - M.getData(): return the pointer to the data space
+  - M(i,j): get an element of the matrix
+
+  The expressions required to be valid for a class M implementing the "mutable matrix concept" are (in addition to those of "matrix concept"):
+  - M = X: assignment operator
+  - M(i,j) = d: assign an element of the matrix
+  - M.setAll(d): set all the elements of the matrix
+*/
+
+//! A full matrix, implements the "mutable matrix concept"
+/*! Owns the data space for the elements */
+class Matrix
+{
+ private:
+  //! Stored in column-major order, as in Fortran and MATLAB
+  double *data;
+  const size_t rows, cols;
+ public:
+  Matrix(size_t rows_arg, size_t cols_arg);
+  Matrix(size_t size_arg);
+  Matrix(const Matrix &arg);
+  ~Matrix();
+  inline size_t getRows() const { return rows; }
+  inline size_t getCols() const { return cols; }
+  inline size_t getLd() const { return rows; }
+  inline double *getData() { return data; }
+  inline const double *getData() const { return data; }
+  inline void setAll(double val) { std::fill_n(data, rows*cols, val); }
+  inline double &operator() (size_t i, size_t j) { return data[i+j*rows]; }
+  inline const double &operator() (size_t i, size_t j) const { return data[i+j*rows]; }
+  //! Assignment operator, only works for matrices of same dimension
+  template<class Mat>
+  Matrix &
+  operator= (const Mat &arg)
+  {
+    assert(rows == arg.getRows() && cols == arg.getCols());
+    for (size_t j = 0; j < cols; j++)
+      memcpy(data + j*rows, arg.getData() + j*arg.getLd(), rows*sizeof(double));
+    return *this;
+  }
+  //! The copy assignment operator, which is not generated by the template assignment operator
+  /*! See C++ standard, §12.8.9, in the footnote */
+  Matrix &operator= (const Matrix &arg);
+};
+
+//! A contiguous submatrix of another matrix, implements the "mutable matrix concept"
+/*! Does not own the data space for the elements, so depends on another matrix */
+class MatrixView
+{
+ private:
+  double *const data;
+  const size_t rows, cols, ld;
+ public:
+  MatrixView(double *data_arg, size_t rows_arg, size_t cols_arg, size_t ld_arg);
+  MatrixView(Matrix &arg, size_t row_offset, size_t col_offset,
+             size_t rows_arg, size_t cols_arg);
+  inline size_t getRows() const { return rows; }
+  inline size_t getCols() const { return cols; }
+  inline size_t getLd() const { return ld; }
+  inline double *getData() { return data; }
+  inline const double *getData() const { return data; }
+  inline void setAll(double val)
+  {
+    for (double *p = data; p < data + rows*ld; p += ld)
+      std::fill_n(p, cols, val);
+  }
+  inline double &operator() (size_t i, size_t j) { return data[i+j*ld]; }
+  inline const double &operator() (size_t i, size_t j) const { return data[i+j*ld]; }
+  //! Assignment operator, only works for matrices of same dimension
+  template<class Mat>
+  MatrixView &
+  operator= (const Mat &arg)
+  {
+    assert(rows == arg.getRows() && cols == arg.getCols());
+    for (size_t j = 0; j < cols; j++)
+      memcpy(data + j*ld, arg.getData() + j*arg.getLd(), rows*sizeof(double));
+    return *this;
+  }
+  //! The copy assignment operator, which is not generated by the template assignment operator
+  /*! See C++ standard, §12.8.9, in the footnote */
+  MatrixView &operator= (const MatrixView &arg);
+};
+
+//! Like MatrixView, but cannot be modified (implements the "matrix concept")
+class MatrixConstView
+{
+ private:
+  const double *const data;
+  const size_t rows, cols, ld;
+ public:
+  MatrixConstView(const double *data_arg, size_t rows_arg, size_t cols_arg, size_t ld_arg);
+  MatrixConstView(const Matrix &arg, size_t row_offset, size_t col_offset,
+                  size_t rows_arg, size_t cols_arg);
+  inline size_t getRows() const { return rows; }
+  inline size_t getCols() const { return cols; }
+  inline size_t getLd() const { return ld; }
+  inline const double *getData() const { return data; }
+  inline const double &operator() (size_t i, size_t j) const { return data[i+j*ld]; }
+};
+
+std::ostream &operator<<(std::ostream &out, const Matrix &M);
+std::ostream &operator<<(std::ostream &out, const MatrixView &M);
+std::ostream &operator<<(std::ostream &out, const MatrixConstView &M);
+
+namespace mat
+{
+  template<class Mat>
+  void
+  print(std::ostream &out, const Mat &M)
+  {
+    for(size_t i = 0; i < M.getRows(); i++)
+      {
+        for(size_t j = 0; j < M.getCols(); j++)
+          out << std::setw(13) << std::right << M(i, j) << " ";
+        out << std::endl;
+      }
+  }
+
+  template<class Mat>
+  inline VectorView
+  get_col(Mat &M, size_t j)
+  {
+    return VectorView(M.getData()+j*M.getLd(), M.getRows(), 1);
+  }
+
+  template<class Mat>
+  inline VectorView
+  get_row(Mat &M, size_t i)
+  {
+    return VectorView(M.getData()+i, M.getCols(), M.getLd());
+  }
+
+  template<class Mat>
+  inline VectorConstView
+  get_col(const Mat &M, size_t j)
+  {
+    return VectorView(M.getData()+j*M.getLd(), M.getRows(), 1);
+  }
+
+  template<class Mat>
+  inline VectorConstView
+  get_row(const Mat &M, size_t i)
+  {
+    return VectorView(M.getData()+i, M.getCols(), M.getLd());
+  }
+
+  template<class Mat1, class Mat2>
+  inline void
+  col_copy(const Mat1 &src, size_t col_src, Mat2 &dest, size_t col_dest)
+  {
+    assert(src.getRows() == dest.getRows()
+           && col_src < src.getCols() && col_dest < dest.getCols());
+    memcpy(dest.getData() + col_dest*dest.getLd(),
+           const_cast<double *>(src.getData()) + col_src*src.getLd(),
+           src.getRows()*sizeof(double));
+  }
+
+  template<class Mat1, class Mat2>
+  inline void
+  col_copy(const Mat1 &src, size_t col_src, size_t row_offset_src, size_t row_nb,
+           Mat2 &dest, size_t col_dest, size_t row_offset_dest)
+  {
+    assert(src.getRows() == dest.getRows()
+           && col_src < src.getCols() && col_dest < dest.getCols()
+           && row_offset_src < src.getRows() && row_offset_src+row_nb <= src.getRows()
+           && row_offset_dest < dest.getRows() && row_offset_dest+row_nb <= dest.getRows());
+    memcpy(dest.getData() + row_offset_dest + col_dest*dest.getLd(),
+           src.getData() + row_offset_src + col_src*src.getLd(),
+           row_nb*sizeof(double));
+  }
+
+  template<class Mat1, class Mat2>
+  inline void
+  row_copy(const Mat1 &src, size_t row_src, Mat2 &dest, size_t row_dest)
+  {
+    assert(src.getCols() == dest.getCols()
+           && row_src < src.getRows() && row_dest < dest.getRows());
+    const double *p1 = src.getData() + row_src;
+    double *p2 = dest.getData() + row_dest;
+    while (p1 < src.getData() + src.getRows() * src.getLd())
+      {
+        *p2 = *p1;
+        p1 += src.getLd();
+        p2 += dest.getLd();
+      }
+  }
+
+  template<class Mat>
+  inline void
+  col_set(Mat &M, size_t col, size_t row_offset, size_t row_nb, double val)
+  {
+    assert(col < M.getCols());
+    assert(row_offset < M.getRows() && row_offset + row_nb <= M.getRows());
+    std::fill_n(M.getData() + M.getLd()*col + row_offset, row_nb, val);
+  }
+
+  //! Copy under the diagonal the elements above the diagonal
+  template<class Mat>
+  inline void
+  copy_upper_to_lower(Mat &M)
+  {
+    size_t d = std::min(M.getCols(), M.getRows());
+    for (size_t i = 0; i < d; i++)
+      for (size_t j = 0; j < i; j++)
+        M(i, j) = M(j, i);
+  }
+
+  //! Copy above the diagonal the elements under the diagonal
+  template<class Mat>
+  inline void
+  copy_lower_to_upper(Mat &M)
+  {
+    size_t d = std::min(M.getCols(), M.getRows());
+    for (size_t i = 0; i < d; i++)
+      for (size_t j = 0; j < i; j++)
+        M(j, i) = M(i, j);
+  }
+
+  //! Fill the matrix with the identity matrix
+  template<class Mat>
+  inline void
+  set_identity(Mat &M)
+  {
+    M.setAll(0.0);
+    size_t d = std::min(M.getCols(), M.getRows());
+    for (size_t i = 0; i < d; i++)
+      M(i, i) = 1.0;
+  }
+
+  //! In-place transpose of a square matrix
+  template<class Mat>
+  inline void
+  transpose(Mat &M)
+  {
+    assert(M.getRows() == M.getCols());
+    for (size_t i = 0; i < M.getRows(); i++)
+      for (size_t j = 0; j < i; j++)
+        std::swap(M(i,j), M(j,i));
+  }
+
+  //! Computes m1 = m1 + m2
+  template<class Mat1, class Mat2>
+  void
+  add(Mat1 m1, const Mat2 m2)
+  {
+    assert(m1.getRows() == m2.getRows() && m1.getCols() == m2.getCols());
+    double *p1 = m1.getData();
+    const double *p2 = m2.getData();
+    while (p1 < m1.getData() + m1.getRows() * m1.getLd())
+      {
+        double *pp1, *pp2;
+        while (pp1 < p1 + m1.getCols())
+          *pp1++ += *pp2++;
+
+        p1 += m1.getLd();
+        p2 += m2.getLd();
+      }
+  }
+
+  //! Computes m1 = m1  m2
+  template<class Mat1, class Mat2>
+  void
+  sub(Mat1 m1, const Mat2 m2)
+  {
+    assert(m1.getRows() == m2.getRows() && m1.getCols() == m2.getCols());
+    double *p1 = m1.getData();
+    const double *p2 = m2.getData();
+    while (p1 < m1.getData() + m1.getRows() * m1.getLd())
+      {
+        double *pp1, *pp2;
+        while (pp1 < p1 + m1.getCols())
+          *pp1++ -= *pp2++;
+
+        p1 += m1.getLd();
+        p2 += m2.getLd();
+      }
+  }
+
+  //! Does m = -m
+  template<class Mat>
+  void
+  negate(Mat m)
+  {
+    double *p = m.getData();
+    while (p < m.getData() + m.getRows() * m.getLd())
+      {
+        double *pp;
+        while (pp < p + m.getCols())
+          {
+            *pp = -*pp;
+            pp++;
+          }
+
+        p += m.getLd();
+      }
+  }
+
+} // End of namespace
+
+#endif
diff --git a/mex/sources/estimation/libmat/QRDecomposition.cc b/mex/sources/estimation/libmat/QRDecomposition.cc
new file mode 100644
index 0000000000000000000000000000000000000000..de7a9466c5f33413e5c073d40ff955cf0ab8b1cb
--- /dev/null
+++ b/mex/sources/estimation/libmat/QRDecomposition.cc
@@ -0,0 +1,36 @@
+/*
+ * Copyright (C) 2010 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 <algorithm> // For std::min
+
+#include "QRDecomposition.hh"
+
+QRDecomposition::QRDecomposition(size_t rows_arg, size_t cols_arg) :
+  rows(rows_arg), cols(cols_arg), mind(std::min(rows, cols)), lwork(rows*cols),
+  H(rows), Q2(rows), v(rows)
+{
+  work = new double[lwork];
+  tau = new double[mind];
+}
+
+QRDecomposition::~QRDecomposition()
+{
+  delete[] work;
+  delete[] tau;
+}
diff --git a/mex/sources/estimation/libmat/QRDecomposition.hh b/mex/sources/estimation/libmat/QRDecomposition.hh
new file mode 100644
index 0000000000000000000000000000000000000000..300a274acb6d41dd22ca4aed9d046232781fd5a5
--- /dev/null
+++ b/mex/sources/estimation/libmat/QRDecomposition.hh
@@ -0,0 +1,70 @@
+/*
+ * Copyright (C) 2010 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 <algorithm> // For std::min()
+
+#include <dynlapack.h>
+
+#include "Vector.hh"
+#include "Matrix.hh"
+#include "BlasBindings.hh"
+
+class QRDecomposition
+{
+private:
+  const size_t rows, cols, mind;
+  lapack_int lwork;
+  double *work, *tau;
+  Matrix H, Q2;
+  Vector v;
+public:
+  /*!
+    \todo Replace heuristic choice for workspace size by a query to determine the optimal size
+   */
+  QRDecomposition(size_t rows_arg, size_t cols_arg);
+  ~QRDecomposition();
+  //! Performs the QR decomposition of a matrix, and multiplies another matrix by Q
+  /*!
+    \param[in,out] A On input, the matrix to be decomposed. On output, equals to the output of dgeqrf
+    \param[in] side Specifies on which side is Q in the multiplication, either "L" or "R"
+    \param[in] trans Specifies whether Q should be transposed before the multiplication, either "T" or "N"
+    \param[in,out] C The matrix to be multiplied by Q, modified in place
+  */
+  template<class Mat1, class Mat2>
+  void computeAndMultByQ(Mat1 &A, const char *side, const char *trans, Mat2 &C);
+};
+
+template<class Mat1, class Mat2>
+void
+QRDecomposition::computeAndMultByQ(Mat1 &A, const char *side, const char *trans, Mat2 &C)
+{
+  assert(A.getRows() == rows && A.getCols() == cols);
+  assert(C.getRows() == rows);
+
+  lapack_int m = rows, n = cols, lda = A.getLd();
+  lapack_int info;
+  dgeqrf(&m, &n, A.getData(), &lda, tau, work, &lwork, &info);
+  assert(info == 0);
+
+  n = C.getCols();
+  lapack_int k = mind, ldc = C.getLd();
+  dormqr(side, trans, &m, &n, &k, A.getData(), &lda, tau, C.getData(), &ldc,
+         work, &lwork, &info);
+  assert(info == 0);
+}
diff --git a/mex/sources/estimation/libmat/Vector.cc b/mex/sources/estimation/libmat/Vector.cc
new file mode 100644
index 0000000000000000000000000000000000000000..f22993e7fc9c90985d97c58c8d26470879612bca
--- /dev/null
+++ b/mex/sources/estimation/libmat/Vector.cc
@@ -0,0 +1,98 @@
+/*
+ * Copyright (C) 2010 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 "Vector.hh"
+
+#include <cstring>
+#include <cassert>
+
+Vector::Vector(size_t size_arg) : size(size_arg)
+{
+  data = new double[size];
+}
+
+Vector::Vector(const Vector &arg) : size(arg.size)
+{
+  memcpy(data, arg.data, size*sizeof(double));
+}
+
+Vector::~Vector()
+{
+  delete[] data;
+}
+
+Vector &
+Vector::operator= (const Vector &arg)
+{
+  assert(size == arg.size);
+  memcpy(data, arg.data, size*sizeof(double));
+  return *this;
+}
+
+std::ostream &
+operator<<(std::ostream &out, const Vector &V)
+{
+  vec::print(out, V);
+  return out;
+}
+
+VectorView::VectorView(Vector &arg, size_t offset, size_t size_arg)
+  : data(arg.getData() + offset), size(size_arg), stride(1)
+{
+}
+
+VectorView::VectorView(double *data_arg, size_t size_arg, size_t stride_arg)
+  : data(data_arg), size(size_arg), stride(stride_arg)
+{
+}
+
+VectorView &
+VectorView::operator= (const VectorView &arg)
+{
+  assert(size == arg.getSize());
+  double *p1;
+  const double *p2;
+  for (p1 = data, p2 = arg.getData(); p1 < data + size*stride; p1 += stride, p2 += arg.getStride())
+    *p1 = *p2;
+  return *this;
+}
+
+std::ostream &
+operator<<(std::ostream &out, const VectorView &V)
+{
+  vec::print(out, V);
+  return out;
+}
+
+VectorConstView::VectorConstView(const Vector &arg, size_t offset, size_t size_arg)
+  : data(arg.getData() + offset), size(size_arg), stride(1)
+{
+}
+
+VectorConstView::VectorConstView(const double *data_arg, size_t size_arg, size_t stride_arg)
+  : data(data_arg), size(size_arg), stride(stride_arg)
+{
+}
+
+std::ostream &
+operator<<(std::ostream &out, const VectorConstView &V)
+{
+  vec::print(out, V);
+  return out;
+}
diff --git a/mex/sources/estimation/libmat/Vector.hh b/mex/sources/estimation/libmat/Vector.hh
new file mode 100644
index 0000000000000000000000000000000000000000..43efcf0bbeda1f7febbca96378bea5b6399dd850
--- /dev/null
+++ b/mex/sources/estimation/libmat/Vector.hh
@@ -0,0 +1,206 @@
+/*
+ * Copyright (C) 2010 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/>.
+ */
+
+#ifndef _VECTOR_HH
+#define _VECTOR_HH
+
+#include <algorithm>
+#include <iostream>
+
+#include <cstdlib>
+#include <cassert>
+
+/*
+  This header defines three vector classes, which implement a "vector concept"
+  (much like the concepts of the Standard Template Library or of Boost
+  Library). The first class is a vector owning the data space for its elements,
+  and the other two are vector "views" of something else (either a subvector of
+  another vector, or a row or column of a matrix). This design philosophy is
+  close to the design of the GNU Scientific Library, but here using the
+  syntactic power of C++ class and templates, while achieving very high
+  efficiency.
+
+  These classes can be used with various templated functions, including
+  wrappers around the BLAS primitives.
+
+  The expressions required to be valid for a class V implementing the "vector concept" are:
+  - V.getSize(): return number of elements
+  - V.getStride(): return the stride, i.e. the offset between two elements in the data space
+  - V.getData(): return the pointer to the data space
+  - V(i): get an element of the vector
+
+  The expressions required to be valid for a class V implementing the "mutable vector concept" are (in addition to those of "vector concept"):
+  - V = X: assignment operator
+  - V(i) = d: assign an element of the vector
+  - V.setAll(d): set all the elements of the vector
+
+*/
+
+//! A full vector, implements the "mutable vector concept"
+/*! Owns the data space for the elements */
+class Vector
+{
+ private:
+  double *data;
+  const size_t size;
+ public:
+  Vector(size_t size_arg);
+  Vector(const Vector &arg);
+  ~Vector();
+  inline size_t getSize() const { return size; }
+  inline size_t getStride() const { return 1; }
+  inline double *getData() { return data; }
+  inline const double *getData() const { return data; }
+  inline void setAll(double val) { std::fill_n(data, size, val); }
+  inline double &operator() (size_t i) { return data[i]; }
+  inline const double &operator() (size_t i) const { return data[i]; }
+  //! Assignment operator, only works for vectors of same size
+  template<class Vec>
+  Vector &operator= (const Vec &arg)
+  {
+    assert(size == arg.getSize());
+    const double *p2 = arg.getData();
+    for(double *p1 = data; p1 < data + size; p1++)
+      {
+        *p1 = *p2;
+        p2 += arg.getStride();
+      }
+    return *this;
+  }
+  //! The copy assignment operator, which is not generated by the template assignment operator
+  /*! See C++ standard, §12.8.9, in the footnote */
+  Vector &operator= (const Vector &arg);
+};
+
+//! A vector view (i.e. a subvector of another vector, or a row or column of a matrix), implements the "mutable vector concept"
+/*! Does not own the data space for the elements, so depends on another vector or matrix */
+class VectorView
+{
+private:
+  double *const data;
+  const size_t size, stride;
+public:
+  VectorView(Vector &arg, size_t offset, size_t size_arg);
+  VectorView(double *data_arg, size_t size_arg, size_t stride_arg);
+  inline size_t getSize() const { return size; }
+  inline size_t getStride() const { return stride; }
+  inline double *getData() { return data; }
+  inline const double *getData() const { return data; }
+  inline void setAll(double val)
+  {
+    for (double *p = data; p < data + size*stride; p += stride)
+      *p = val;
+  }
+  inline double &operator() (size_t i) { return data[i*stride]; }
+  inline const double &operator() (size_t i) const { return data[i*stride]; }
+  //! Assignment operator, only works for vectors of same size
+  template<class Vec>
+  VectorView &
+  operator= (const Vec &arg)
+  {
+    assert(size == arg.getSize());
+    double *p1;
+    const double *p2;
+    for (p1 = data, p2 = arg.getData(); p1 < data + size*stride; p1 += stride, p2 += arg.getStride())
+      *p1 = *p2;
+    return *this;
+  }
+  //! The copy assignment operator, which is not generated by the template assignment operator
+  /*! See C++ standard, §12.8.9, in the footnote */
+  VectorView &operator= (const VectorView &arg);
+};
+
+//! Like a VectorView, but the data cannot be modified (implements the "vector concept")
+class VectorConstView
+{
+private:
+  const double *const data;
+  const size_t size, stride;
+public:
+  VectorConstView(const Vector &arg, size_t offset, size_t size_arg);
+  VectorConstView(const double *data_arg, size_t size_arg, size_t stride_arg);
+  inline size_t getSize() const { return size; }
+  inline size_t getStride() const { return stride; }
+  inline const double *getData() const { return data; }
+  inline const double &operator() (size_t i) const { return data[i*stride]; }
+};
+
+std::ostream &operator<<(std::ostream &out, const Vector &V);
+std::ostream &operator<<(std::ostream &out, const VectorView &V);
+std::ostream &operator<<(std::ostream &out, const VectorConstView &V);
+
+namespace vec
+{
+  template<class Vec>
+  void
+  print(std::ostream &out, const Vec &V)
+  {
+    for (size_t i = 0; i < V.getSize(); i++)
+      out << V(i) << " ";
+    out << std::endl;
+  }
+
+  //! Computes v1 = v1 + v2
+  template<class Vec1, class Vec2>
+  void
+  add(Vec1 v1, const Vec2 v2)
+  {
+    assert(v1.getSize() == v2.getSize());
+    double *p1 = v1.getData();
+    const double *p2 = v2.getData();
+    while (p1 < v1.getData() + v1.getSize() * v1.getStride())
+      {
+        *p1 += *p2;
+        p1 += v1.getStride();
+        p2 += v2.getStride();
+      }
+  }
+
+  //! Computes v1 = v1 - v2
+  template<class Vec1, class Vec2>
+  void
+  sub(Vec1 v1, const Vec2 v2)
+  {
+    assert(v1.getSize() == v2.getSize());
+    double *p1 = v1.getData();
+    const double *p2 = v2.getData();
+    while (p1 < v1.getData() + v1.getSize() * v1.getStride())
+      {
+        *p1 -= *p2;
+        p1 += v1.getStride();
+        p2 += v2.getStride();
+      }
+  }
+
+  //! Does v = -v
+  template<class Vec>
+  void
+  negate(Vec v)
+  {
+    double *p = v.getData();
+    while (p < v.getData() + v.getSize() * v.getStride())
+      {
+        *p = -*p;
+        p += v.getStride();
+      }
+  }
+
+} // End of namespace
+
+#endif
diff --git a/mex/sources/estimation/libmat/tests/Makefile.am b/mex/sources/estimation/libmat/tests/Makefile.am
new file mode 100644
index 0000000000000000000000000000000000000000..353cd327fc459c7f653771a58a27b7cb97d14491
--- /dev/null
+++ b/mex/sources/estimation/libmat/tests/Makefile.am
@@ -0,0 +1,9 @@
+check_PROGRAMS = test-qr test-gsd
+
+test_qr_SOURCES = ../Matrix.cc ../Vector.cc ../QRDecomposition.cc test-qr.cc
+test_qr_LDADD = $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS)
+test_qr_CPPFLAGS = -I.. -I../../../
+
+test_gsd_SOURCES = ../Matrix.cc ../Vector.cc ../GeneralizedSchurDecomposition.cc test-gsd.cc
+test_gsd_LDADD = $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS)
+test_gsd_CPPFLAGS = -I.. -I../../../
diff --git a/mex/sources/estimation/libmat/tests/test-gsd.cc b/mex/sources/estimation/libmat/tests/test-gsd.cc
new file mode 100644
index 0000000000000000000000000000000000000000..52211076884205fe9fe2be821a6afdfc885f1fe9
--- /dev/null
+++ b/mex/sources/estimation/libmat/tests/test-gsd.cc
@@ -0,0 +1,61 @@
+/*
+ * Copyright (C) 2010 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 <iostream>
+
+#include "GeneralizedSchurDecomposition.hh"
+
+int
+main(int argc, char **argv)
+{
+  size_t n = 3;
+  double D_data[] = { 1, 2, 3,
+                      4, 5, 6,
+                      7, 8, 9 };
+  double E_data[] = { 1, -3, 4,
+                      -7, 9, 1,
+                      -3, 4, 0 };
+  MatrixView D(D_data, n, n, n), E(E_data, n, n, n);
+
+  // Need to transpose because internally matrices are in column-major order
+  mat::transpose(D);
+  mat::transpose(E);
+
+  std::cout << "Matrix D: " << std::endl << D << std::endl;
+  std::cout << "Matrix E: " << std::endl << E << std::endl;
+
+  GeneralizedSchurDecomposition GSD(n, 1.00001);
+
+  Matrix S(n), T(n), Z(n);
+  size_t sdim;
+
+  GSD.compute(D, E, S, T, Z, sdim);
+
+  std::cout << "Matrix S: " << std::endl << S << std::endl;
+  std::cout << "Matrix T: " << std::endl << T << std::endl;
+  std::cout << "Matrix Z: " << std::endl << Z << std::endl;
+
+  Vector eig_real(n), eig_cmplx(n);
+  GSD.getGeneralizedEigenvalues(eig_real, eig_cmplx);
+
+  std::cout << "Real part of generalized eigenvalues: " << std::endl << eig_real << std::endl;
+  std::cout << "Complex part of generalized eigenvalues: " << std::endl << eig_cmplx << std::endl;
+
+  return 0;
+}
diff --git a/mex/sources/estimation/libmat/tests/test-qr.cc b/mex/sources/estimation/libmat/tests/test-qr.cc
new file mode 100644
index 0000000000000000000000000000000000000000..ad25e142a6adebb0a9d96d1371aef289d6eb04a5
--- /dev/null
+++ b/mex/sources/estimation/libmat/tests/test-qr.cc
@@ -0,0 +1,56 @@
+/*
+ * Copyright (C) 2010 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 <iostream>
+
+#include "BlasBindings.hh"
+#include "QRDecomposition.hh"
+
+int
+main(int argc, char **argv)
+{
+  size_t m = 4, n = 3;
+  Matrix S(m, n), Q(m), A(m, n), B(m);
+  QRDecomposition QRD(m, n);
+
+  for (size_t i = 0; i < m; i++)
+    for (size_t j = 0; j < n; j++)
+      S(i, j) = i*n + j + 1;
+
+  std::cout << "Matrix to be decomposed:" << std::endl << S << std::endl;
+
+  mat::set_identity(Q);
+
+  QRD.computeAndMultByQ(S, "L", "N", Q);
+
+  std::cout << "Matrix Q:" << std::endl << Q << std::endl;
+
+  blas::gemm("T", "N", 1.0, Q, Q, 0.0, B);
+
+  std::cout << "Matrix Q'*Q:" << std::endl << B << std::endl;
+
+  for (size_t j = 0; j < n; j++)
+    mat::col_set(S, j, j+1, m-j-1, 0);
+
+  std::cout << "Matrix R:" << std::endl << S << std::endl;
+
+  blas::gemm("N", "N", 1.0, Q, S, 0.0, A);
+
+  std::cout << "Product Q*R:" << std::endl << A << std::endl;
+}
diff --git a/mex/sources/estimation/tests/Makefile.am b/mex/sources/estimation/tests/Makefile.am
new file mode 100644
index 0000000000000000000000000000000000000000..2099edfb7ddd568223480f25e3748fd6ea0e55e3
--- /dev/null
+++ b/mex/sources/estimation/tests/Makefile.am
@@ -0,0 +1,5 @@
+check_PROGRAMS = test-dr
+
+test_dr_SOURCES = ../libmat/Matrix.cc ../libmat/Vector.cc ../libmat/QRDecomposition.cc ../libmat/GeneralizedSchurDecomposition.cc ../libmat/LUSolver.cc ../DecisionRules.cc test-dr.cc
+test_dr_LDADD = $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS)
+test_dr_CPPFLAGS = -I.. -I../libmat -I../../
diff --git a/mex/sources/estimation/tests/test-dr.cc b/mex/sources/estimation/tests/test-dr.cc
new file mode 100644
index 0000000000000000000000000000000000000000..373794cf578dedfb4871fce788d51632234a331d
--- /dev/null
+++ b/mex/sources/estimation/tests/test-dr.cc
@@ -0,0 +1,25 @@
+/*
+ * Copyright (C) 2010 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 "DecisionRules.hh"
+
+int
+main(int argc, char **argv)
+{
+}