mex/sources/estimation/libmat/Matrix.cc

-/*
- * 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 "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)
-{
-}
-
-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)
-{
-}
-
-std::ostream &
-operator<<(std::ostream &out, const MatrixConstView &M)
-{
-  mat::print(out, M);
-  return out;
-}

mex/sources/estimation/libmat/Matrix.hh

-/*
- * Copyright (C) 2010-2012 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 <cmath>
-#include <vector>
-
-#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);
-  virtual ~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);
-  template<class Mat>
-  MatrixView(Mat &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());
-  }
-  virtual ~MatrixView(){};
-  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 + cols*ld; p += ld)
-      std::fill_n(p, rows, 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);
-  template<class Mat>
-  MatrixConstView(const Mat &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());
-  }
-  virtual ~MatrixConstView(){};
-  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
-{
-  //define nullVec (const vector<int>(0)) for assign and order by vector
-  // It is used as a proxy for the ":" matlab operator:
-  // i.e. zero sized int vector, nullVec, is interpreted as if one supplied ":"
-  const std::vector<size_t> nullVec(0);
-
-  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 VectorConstView(M.getData()+j*M.getLd(), M.getRows(), 1);
-  }
-
-  template<class Mat>
-  inline VectorConstView
-  get_row(const Mat &M, size_t i)
-  {
-    return VectorConstView(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(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.getCols() * 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 = M2' (even for rectangular matrices)
-  template<class Mat1, class Mat2>
-  inline void
-  transpose(Mat1 &M1, const Mat2 &M2)
-  {
-    assert(M1.getRows() == M2.getCols() && M1.getCols() == M2.getRows());
-    for (size_t i = 0; i < M1.getRows(); i++)
-      for (size_t j = 0; j < M1.getCols(); j++)
-        M1(i, j) = M2(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.getCols() * m1.getLd())
-      {
-        double *pp1 = p1;
-        const double *pp2 = p2;
-        while (pp1 < p1 + m1.getRows())
-          *pp1++ += *pp2++;
-
-        p1 += m1.getLd();
-        p2 += m2.getLd();
-      }
-  }
-
-  //! Computes m1 = m1 + number
-  template<class Mat1>
-  void
-  add(Mat1 &m1, double d)
-  {
-    double *p1 = m1.getData();
-
-    while (p1 < m1.getData() + m1.getCols() * m1.getLd())
-      {
-        double *pp1 = p1;
-        while (pp1 < p1 + m1.getRows())
-          *pp1++ += d;
-
-        p1 += m1.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.getCols() * m1.getLd())
-      {
-        double *pp1 = p1;
-        const double *pp2 = p2;
-        while (pp1 < p1 + m1.getRows())
-          *pp1++ -= *pp2++;
-
-        p1 += m1.getLd();
-        p2 += m2.getLd();
-      }
-  }
-
-  //! Computes m1 = m1 - number
-  template<class Mat1>
-  void
-  sub(Mat1 &m1, double d)
-  {
-    add(m1, -1.0*d);
-  }
-
-  //! Does m = -m
-  template<class Mat>
-  void
-  negate(Mat &m)
-  {
-    double *p = m.getData();
-    while (p < m.getData() + m.getCols() * m.getLd())
-      {
-        double *pp = p;
-        while (pp < p + m.getRows())
-          {
-            *pp = -*pp;
-            pp++;
-          }
-
-        p += m.getLd();
-      }
-  }
-
-  // Computes the infinite norm of a matrix
-  template<class Mat>
-  double
-  nrminf(const Mat &m)
-  {
-    double nrm = 0;
-    const double *p = m.getData();
-    while (p < m.getData() + m.getCols() * m.getLd())
-      {
-        const double *pp = p;
-        while (pp < p + m.getRows())
-          {
-            if (fabs(*pp) > nrm)
-              nrm = fabs(*pp);
-            pp++;
-          }
-
-        p += m.getLd();
-      }
-    return nrm;
-  }
-
-  // emulates Matlab command A(:,b)=B(:,d) where b,d are size_t vectors or nullVec as a proxy for ":")
-  // i.e. zero sized vector (or mat::nullVec) is interpreted as if one supplied ":" in matlab
-  template<class Mat1, class Mat2>
-  void
-  reorderColumnsByVectors(Mat1 &a, const std::vector<size_t> &vToCols,
-                          const Mat2 &b, const std::vector<size_t> &vcols)
-  {
-    size_t ncols = 0, toncols = 0;
-    const std::vector<size_t> *vpToCols = 0, *vpCols = 0;
-    std::vector<size_t> tmpvpToCols(0), tmpvpCols(0);
-    assert(b.getRows() ==  a.getRows());
-
-    if (vToCols.size() == 0  && vcols.size() == 0)
-      a = b;
-    else
-      {
-        if (vToCols.size() == 0)
-          {
-            toncols = a.getCols();
-            tmpvpToCols.reserve(toncols);
-            for (size_t i = 0; i < toncols; ++i)
-              tmpvpToCols[i] = i;
-            vpToCols = (const std::vector<size_t> *)&tmpvpToCols;
-          }
-        else
-          {
-            for (size_t i = 0; i < vToCols.size(); ++i)
-              {
-                assert(vToCols[i] < a.getCols()); //Negative or too large indices
-                toncols++;
-              }
-            assert(toncols <= a.getCols()); // check wrong dimensions for assignment by vector
-            vpToCols = &vToCols;
-          }
-
-        if (vcols.size() == 0)
-          {
-            ncols = b.getCols();
-            tmpvpCols.reserve(ncols);
-            for (size_t i = 0; i < ncols; ++i)
-              tmpvpCols[i] = i;
-            vpCols = (const std::vector<size_t> *)&tmpvpCols;
-          }
-        else
-          {
-            for (size_t i = 0; i < vcols.size(); ++i)
-              {
-                assert(vcols[i] < b.getCols()); //Negative or too large indices
-                ncols++;
-              }
-            assert(ncols <= b.getCols()); // check wrong dimensions for assignment by vector
-            vpCols = &vcols;
-          }
-
-        assert(toncols == ncols && ncols > 0);
-        for (size_t j = 0; j < ncols; ++j)
-          col_copy(b,  (*vpCols)[j], a,  (*vpToCols)[j]);
-      }
-  }
-
-  // emulates Matlab command A(a,:)=B(c,:) where a,c are size_t vectors or nullVec as a proxy for ":")
-  // i.e. zero sized vector (or mat::nullVec) is interpreted as if one supplied ":" in matlab
-  template<class Mat1, class Mat2>
-  void
-  reorderRowsByVectors(Mat1 &a, const std::vector<size_t> &vToRows,
-                       const Mat2 &b, const std::vector<size_t> &vrows)
-  {
-    size_t nrows = 0, tonrows = 0;
-    const std::vector<size_t>  *vpToRows = 0, *vpRows = 0;
-    std::vector<size_t>  tmpvpToRows(0), tmpvpRows(0);
-
-    //assert(b.getRows() >=  a.getRows() && b.getCols() ==  a.getCols());
-    assert(b.getCols() ==  a.getCols());
-    if (vToRows.size() == 0  && vrows.size() == 0)
-      a = b;
-    else
-      {
-        if (vToRows.size() == 0)
-          {
-            tonrows = a.getRows();
-            tmpvpToRows.reserve(tonrows);
-            for (size_t i = 0; i < tonrows; ++i)
-              tmpvpToRows[i] = i;
-            vpToRows = (const std::vector<size_t> *)&tmpvpToRows;
-          }
-        else
-          {
-            for (size_t i = 0; i < vToRows.size(); ++i)
-              {
-                assert(vToRows[i] < a.getRows()); //Negative or too large indices
-                tonrows++;
-              }
-            assert(tonrows <= a.getRows()); // check wrong dimensions for assignment by vector
-            vpToRows = &vToRows;
-          }
-
-        if (vrows.size() == 0)
-          {
-            nrows = b.getRows();
-            tmpvpRows.reserve(nrows);
-            for (size_t i = 0; i < nrows; ++i)
-              tmpvpRows[i] = i;
-            vpRows = (const std::vector<size_t> *)&tmpvpRows;
-          }
-        else
-          {
-            for (size_t i = 0; i < vrows.size(); ++i)
-              {
-                assert(vrows[i] < b.getRows()); //Negative or too large indices
-                nrows++;
-              }
-            assert(nrows <= b.getRows()); // check wrong dimensions for assignment by vector
-            vpRows = &vrows;
-          }
-
-        assert(tonrows == nrows && nrows > 0);
-        for (size_t i = 0; i < nrows; ++i)
-          row_copy(b, (*vpRows)[i], a, (*vpToRows)[i]);
-      }
-  }
-
-  // emulates Matlab command A(a,b)=B(c,d) where a,b,c,d are size_t vectors or nullVec as a proxy for ":")
-  // i.e. zero sized vector (or mat::nullVec) is interpreted as if one supplied ":" in matlab
-  template<class Mat1, class Mat2>
-  void
-  assignByVectors(Mat1 &a, const std::vector<size_t> &vToRows, const std::vector<size_t> &vToCols,
-                  const Mat2 &b, const std::vector<size_t> &vrows, const std::vector<size_t> &vcols)
-  {
-    size_t nrows = 0, ncols = 0, tonrows = 0, toncols = 0;
-    const std::vector<size_t> *vpToCols = 0, *vpToRows = 0, *vpRows = 0, *vpCols = 0;
-    std::vector<size_t> tmpvpToCols(0), tmpvpToRows(0), tmpvpRows(0), tmpvpCols(0);
-
-    if (vToRows.size() == 0 && vToCols.size() == 0 && vrows.size() == 0 && vcols.size() == 0)
-      a = b;
-    else if (vToRows.size() == 0 && vrows.size() == 0)                                                                                                                                                                                                                                                                                                                                            // just reorder columns
-      reorderColumnsByVectors(a, vToCols, b, vcols);
-    else if (vToCols.size() == 0 && vcols.size() == 0)                                                                                                                                                                                                                                                                                                                                            // just reorder rows
-      reorderRowsByVectors(a, vToRows, b, vrows);
-    else
-      {
-        if (vToRows.size() == 0)
-          {
-            tonrows = a.getRows();
-            tmpvpToRows.reserve(tonrows);
-            for (size_t i = 0; i < tonrows; ++i)
-              tmpvpToRows[i] = i;
-            vpToRows = (const std::vector<size_t> *)&tmpvpToRows;
-          }
-        else
-          {
-            for (size_t i = 0; i < vToRows.size(); ++i)
-              {
-                assert(vToRows[i] < a.getRows()); //Negative or too large indices
-                tonrows++;
-              }
-            assert(tonrows <= a.getRows()); // check wrong dimensions for assignment by vector
-            vpToRows = &vToRows;
-          }
-
-        if (vToCols.size() == 0)
-          {
-            toncols = a.getCols();
-            tmpvpToCols.reserve(toncols);
-            for (size_t i = 0; i < toncols; ++i)
-              tmpvpToCols[i] = i;
-            vpToCols = (const std::vector<size_t> *)&tmpvpToCols;
-          }
-        else
-          {
-            for (size_t i = 0; i < vToCols.size(); ++i)
-              {
-                assert(vToCols[i] < a.getCols()); //Negative or too large indices
-                toncols++;
-              }
-            assert(toncols <= a.getCols()); // check wrong dimensions for assignment by vector
-            vpToCols = &vToCols;
-          }
-
-        if (vrows.size() == 0)
-          {
-            nrows = b.getRows();
-            tmpvpRows.reserve(nrows);
-            for (size_t i = 0; i < nrows; ++i)
-              tmpvpRows[i] = i;
-            vpRows = (const std::vector<size_t> *)&tmpvpRows;
-          }
-        else
-          {
-            for (size_t i = 0; i < vrows.size(); ++i)
-              {
-                assert(vrows[i] < b.getRows()); //Negative or too large indices
-                nrows++;
-              }
-            assert(nrows <= b.getRows()); // check wrong dimensions for assignment by vector
-            vpRows = &vrows;
-          }
-
-        if (vcols.size() == 0)
-          {
-            ncols = b.getCols();
-            tmpvpCols.reserve(ncols);
-            for (size_t i = 0; i < ncols; ++i)
-              tmpvpCols[i] = i;
-            vpCols = (const std::vector<size_t> *)&tmpvpCols;
-          }
-        else
-          {
-            for (size_t i = 0; i < vcols.size(); ++i)
-              {
-                assert(vcols[i] < b.getCols()); //Negative or too large indices
-                ncols++;
-              }
-            assert(ncols <= b.getCols()); // check wrong dimensions for assignment by vector
-            vpCols = &vcols;
-          }
-
-        assert(tonrows == nrows && toncols == ncols && nrows * ncols > 0);
-        for (size_t i = 0; i < nrows; ++i)
-          for (size_t j = 0; j < ncols; ++j)
-            a((*vpToRows)[i], (*vpToCols)[j]) = b((*vpRows)[i], (*vpCols)[j]);
-      }
-  }
-
-  //emulates Matlab repmat: Mat2 = multv*multh tiled [Mat1]
-  template<class Mat1, class Mat2 >
-  void
-  repmat(Mat1 &a, size_t multv, size_t multh, Mat2 &repMat)   // vertical and horisontal replicators
-  {
-    assert(repMat.getRows() == multv * a.getRows() && repMat.getCols() == multh * a.getCols());
-    for (size_t i = 0; i < multv; ++i)
-      for (size_t j = 0; j < multh; ++j)
-        for (size_t k = 0; k < a.getCols(); ++k)
-          col_copy(a, k, 0, a.getRows(), repMat, a.getCols() * j + k, a.getRows() * i);
-  };
-
-  template<class Mat1, class Mat2>
-  bool
-  isDiff(const Mat1 &m1, const Mat2 &m2, const double tol = 0.0)
-  {
-    assert(m2.getRows() == m1.getRows() && m2.getCols() == m1.getCols());
-    const double *p1 = m1.getData();
-    const double *p2 = m2.getData();
-    while (p1 < m1.getData() + m1.getCols() * m1.getLd())
-      {
-        const double *pp1 = p1;
-        const double *pp2 = p2;
-        while (pp1 < p1 + m1.getRows())
-          if (fabs(*pp1++ - *pp2++) > tol)
-            return true;
-
-        p1 += m1.getLd();
-        p2 += m2.getLd();
-      }
-    return false;
-  }
-
-  //traverse the upper triangle only along diagonals where higher changes occur
-  template<class Mat1, class Mat2>
-  bool
-  isDiffSym(const Mat1 &m1, const Mat2 &m2, const double tol = 0.0)
-  {
-    assert(m2.getRows() == m1.getRows() && m2.getCols() == m1.getCols()
-           && m2.getRows() == m1.getCols() && m2.getCols() == m1.getRows());
-    for (size_t i = 0; i < m1.getCols(); i++)
-      for (size_t j = 0; i + j < m1.getCols(); j++)
-        if (fabs(m1(j, j + i) - m2(j, j + i)) > tol)
-          return true;
-    return false;
-  }
-
-} // End of namespace
-
-#endif

mex/sources/estimation/libmat/QRDecomposition.cc

-/*
- * 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, size_t cols2_arg) :
-  rows(rows_arg), cols(cols_arg), mind(std::min(rows, cols)), cols2(cols2_arg),
-  lwork(rows*cols), lwork2(cols2)
-{
-  work = new double[lwork];
-  work2 = new double[lwork2];
-
-  tau = new double[mind];
-}
-
-QRDecomposition::~QRDecomposition()
-{
-  delete[] work;
-  delete[] work2;
-  delete[] tau;
-}

mex/sources/estimation/libmat/QRDecomposition.hh

-/*
- * 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 <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;
-  //! Number of columns of the matrix to be left-multiplied by Q
-  const size_t cols2;
-  lapack_int lwork, lwork2;
-  double *work, *work2, *tau;
-public:
-  /*!
-    \todo Replace heuristic choice for workspace size by a query to determine the optimal size
-    \param[in] rows_arg Number of rows of the matrix to decompose
-    \param[in] cols_arg Number of columns of the matrix to decompose
-    \param[in] cols2_arg Number of columns of the matrix to be multiplied by Q
-  */
-  QRDecomposition(size_t rows_arg, size_t cols_arg, size_t cols2_arg);
-  virtual ~QRDecomposition();
-  //! Performs the QR decomposition of a matrix, and left-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] trans Specifies whether Q should be transposed before the multiplication, either "T" or "N"
-    \param[in,out] C The matrix to be left-multiplied by Q, modified in place
-  */
-  template<class Mat1, class Mat2>
-  void computeAndLeftMultByQ(Mat1 &A, const char *trans, Mat2 &C);
-};
-
-template<class Mat1, class Mat2>
-void
-QRDecomposition::computeAndLeftMultByQ(Mat1 &A, const char *trans, Mat2 &C)
-{
-  assert(A.getRows() == rows && A.getCols() == cols);
-  assert(C.getRows() == rows && C.getCols() == cols2);
-
-  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 = cols2;
-  lapack_int k = mind, ldc = C.getLd();
-  dormqr("L", trans, &m, &n, &k, A.getData(), &lda, tau, C.getData(), &ldc,
-         work2, &lwork2, &info);
-  assert(info == 0);
-}

mex/sources/estimation/libmat/VDVEigDecomposition.cc

-/*
- * 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 "VDVEigDecomposition.hh"
-
-VDVEigDecomposition::VDVEigDecomposition(const Matrix &m) throw (VDVEigException) :
-  lda(m.getLd()), n(m.getCols()), lwork(-1),
-  info(0),   converged(false), V(m), D(n)
-{
-  if (m.getRows() != m.getCols())
-    throw (VDVEigException(info, "Matrix is not square in VDVEigDecomposition constructor"));
-
-  double tmpwork;
-  dsyev("V", "U", &n, V.getData(), &lda, D.getData(), &tmpwork, &lwork, &info);
-  lwork = (lapack_int) tmpwork;
-  work = new double[lwork];
-  if (info < 0)
-    throw (VDVEigException(info, "Internal error in VDVEigDecomposition constructor"));
-}
-
-VDVEigDecomposition::VDVEigDecomposition(size_t inn) throw (VDVEigException) :
-  lda(inn), n(inn), lwork(3*inn-1),
-  info(0),   converged(false), V(inn), D(inn)
-{
-  work = new double[lwork];
-};
-
-std::ostream &
-operator<<(std::ostream &out, const VDVEigDecomposition::VDVEigException &e)
-{
-  out << " info " << e.info << ": " << e.message  << "\n";
-  return out;
-}

mex/sources/estimation/libmat/VDVEigDecomposition.hh

-/*
- * 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/>.
- */
-
-/**
- *  Calculates vector of Eigen values D and matrix of Eigen vectors V
- *  of input matrix m
- */
-
-#if !defined(VDVEigDec_INCLUDE)
-#define VDVEigDec_INCLUDE
-
-#include <cstdlib>
-#include "Vector.hh"
-#include "Matrix.hh"
-#include <dynlapack.h>
-
-class VDVEigDecomposition
-{
-  lapack_int lda, n;
-  lapack_int lwork, info;
-  double *work;
-  bool converged;
-  Matrix V;
-  Vector D;
-public:
-
-public:
-  class VDVEigException
-  {
-  public:
-    const lapack_int info;
-    std::string message;
-    VDVEigException(lapack_int info_arg, std::string message_arg) :
-      info(info_arg), message(message_arg)
-    {
-    };
-  };
-
-  /**
-   *  This constructor only creates optimal workspace using
-   *  input matrix m
-   */
-  VDVEigDecomposition(const Matrix &m) throw (VDVEigException);
-
-  /**
-   *  This constructoro only crates workspace using the size of
-   *  the input matrix m
-   */
-  VDVEigDecomposition(size_t n) throw (VDVEigException);
-
-  virtual ~VDVEigDecomposition()
-  {
-    delete[] work;
-  };
-  template <class Mat>
-  void calculate(const Mat &H) throw (VDVEigException);
-  // get eigenvalues
-  Vector &
-  getD()
-  {
-    return D;
-  };
-  // get eigen vector
-  Matrix &
-  getV()
-  {
-    return V;
-  };
-  // check if converged
-  bool
-  hasConverged()
-  {
-    return converged;
-  };
-};
-
-template <class Mat>
-void
-VDVEigDecomposition::calculate(const Mat &m) throw (VDVEigException)
-{
-  info = 0;
-  if (m.getRows() != m.getCols())
-    throw (VDVEigException(info, "Matrix is not square in VDVEigDecomposition calculate"));
-
-  if (m.getCols() != (size_t) n  || m.getLd() != (size_t) lda)
-    throw (VDVEigException(info, "Matrix not matching VDVEigDecomposition class"));
-
-  double tmpwork;
-  lapack_int tmplwork = -1;
-  V = m;
-  dsyev("V", "U", &n, V.getData(), &lda, D.getData(), &tmpwork, &tmplwork, &info);
-  if (lwork < (lapack_int) tmpwork)
-    {
-      lwork = (lapack_int) tmpwork;
-      delete[] work;
-      work = new double[lwork];
-    }
-  dsyev("V", "U", &n, V.getData(), &lda, D.getData(), work, &lwork, &info);
-
-  if (info < 0)
-    throw (VDVEigException(info, "Internal error in VDVEigDecomposition calculation"));
-  converged = true;
-  if (info)
-    converged = false;
-}
-
-std::ostream &operator<<(std::ostream &out, const VDVEigDecomposition::VDVEigException &e);
-
-#endif

mex/sources/estimation/libmat/Vector.cc

-/*
- * 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 "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)
-{
-  data = new double[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(double *data_arg, size_t size_arg, size_t stride_arg) :
-  data(data_arg), size(size_arg), stride(stride_arg)
-{
-}
-
-VectorView::VectorView(Vector &arg, size_t offset, size_t size_arg) :
-  data(arg.getData() + offset*arg.getStride()), size(size_arg), stride(arg.getStride())
-{
-  assert(offset < arg.getSize() && offset + size <= arg.getSize());
-}
-
-VectorView::VectorView(VectorView &arg, size_t offset, size_t size_arg) :
-  data(arg.getData() + offset*arg.getStride()), size(size_arg), stride(arg.getStride())
-{
-  assert(offset < arg.getSize() && offset + size <= arg.getSize());
-}
-
-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;
-}
-
-VectorConstView::VectorConstView(const double *data_arg, size_t size_arg, size_t stride_arg) :
-  data(data_arg), size(size_arg), stride(stride_arg)
-{
-}
-
-VectorConstView::VectorConstView(const Vector &arg, size_t offset, size_t size_arg) :
-  data(arg.getData() + offset*arg.getStride()), size(size_arg), stride(arg.getStride())
-{
-  assert(offset < arg.getSize() && offset + size <= arg.getSize());
-}
-
-VectorConstView::VectorConstView(const VectorView &arg, size_t offset, size_t size_arg) :
-  data(arg.getData() + offset*arg.getStride()), size(size_arg), stride(arg.getStride())
-{
-  assert(offset < arg.getSize() && offset + size <= arg.getSize());
-}
-
-VectorConstView::VectorConstView(const VectorConstView &arg, size_t offset, size_t size_arg) :
-  data(arg.getData() + offset*arg.getStride()), size(size_arg), stride(arg.getStride())
-{
-  assert(offset < arg.getSize() && offset + size <= arg.getSize());
-}
-
-std::ostream &
-operator<<(std::ostream &out, const VectorView &V)
-{
-  vec::print(out, V);
-  return out;
-}
-
-std::ostream &
-operator<<(std::ostream &out, const VectorConstView &V)
-{
-  vec::print(out, V);
-  return out;
-}

mex/sources/estimation/libmat/Vector.hh

-/*
- * 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/>.
- */
-
-#ifndef _VECTOR_HH
-#define _VECTOR_HH
-
-#include <algorithm>
-#include <iostream>
-
-#include <cstdlib>
-#include <cassert>
-#include <cmath>
-
-/*
-  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);
-  virtual ~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(double *data_arg, size_t size_arg, size_t stride_arg);
-  /* Can't use a template for the 2 constructors below: this would override the
-  constructor which uses a pointer, because the argument list is the same */
-  VectorView(Vector &arg, size_t offset, size_t size_arg);
-  VectorView(VectorView &arg, size_t offset, size_t size_arg);
-  virtual ~VectorView(){};
-  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 double *data_arg, size_t size_arg, size_t stride_arg);
-  /* Can't use a template for the 3 constructors below: this would override the
-     constructor which uses a pointer, because the argument list is the same */
-  VectorConstView(const Vector &arg, size_t offset, size_t size_arg);
-  VectorConstView(const VectorView &arg, size_t offset, size_t size_arg);
-  VectorConstView(const VectorConstView &arg, size_t offset, size_t size_arg);
-
-  virtual ~VectorConstView() {};
-  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();
-      }
-  }
-
-  // Computes the infinite norm of a vector
-  template<class Vec>
-  double
-  nrminf(const Vec &v)
-  {
-    double nrm = 0;
-    const double *p = v.getData();
-    while (p < v.getData() + v.getSize() * v.getStride())
-      {
-        if (fabs(*p) > nrm)
-          nrm = fabs(*p);
-        p += v.getStride();
-      }
-    return nrm;
-  }
-
-  // Computes the sum, min and max of a vector
-  // and returns double mean=sum/size
-  template<class Vec>
-  double
-  meanSumMinMax(double &sum, double &min, double &max, const Vec &v)
-  {
-    sum = 0;
-    min = max = v(0);
-    const double *p = v.getData();
-    while (p < v.getData() + v.getSize() * v.getStride())
-      {
-        if ((*p) > max)
-          max = (*p);
-        if ((*p) < min)
-          min = (*p);
-        sum += *p;
-        p += v.getStride();
-      }
-    return sum/v.getSize();
-  }
-} // End of namespace
-
-#endif

mex/sources/estimation/libmat/tests/Makefile.am

-check_PROGRAMS = test-qr test-gsd test-lu test-repmat
-
-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../../../
-
-test_lu_SOURCES = ../Matrix.cc ../Vector.cc ../LUSolver.cc test-lu.cc
-test_lu_LDADD = $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS)
-test_lu_CPPFLAGS = -I.. -I../../../
-
-test_repmat_SOURCES = ../Matrix.cc ../Vector.cc test-repmat.cc
-test_repmat_CPPFLAGS = -I..
-
-check-local:
-	./test-qr
-	./test-gsd
-	./test-lu
-	./test-repmat

mex/sources/estimation/libmat/tests/test-gsd.cc

-/*
- * 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 << "D =" << std::endl << D << std::endl;
-  std::cout << "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 << "S =" << std::endl << S << std::endl;
-  std::cout << "T =" << std::endl << T << std::endl;
-  std::cout << "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;
-}

mex/sources/estimation/libmat/tests/test-lu.cc

-/*
- * 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 "LUSolver.hh"
-
-int
-main(int argc, char **argv)
-{
-  size_t m = 4, n = 3;
-
-  double A_data[] = { -1, 2, 3,
-                      4, -5, 6,
-                      7, 8, -9 };
-  double B_data[] = { 1, -3, 4, 5,
-                      -7, 9, 1, 7,
-                      -3, 4, 0, -2 };
-  MatrixView A(A_data, n, n, n), B_prime(B_data, m, n, m);
-  Matrix B(n, m);
-
-  mat::transpose(A);
-  mat::transpose(B, B_prime);
-
-  std::cout << "A =" << std::endl << A << std::endl
-            << "B =" << std::endl << B << std::endl;
-
-  LUSolver LU(n);
-
-  LU.invMult("N", A, B);
-
-  std::cout << "A\\B =" << std::endl << B;
-
-  // Check this is the right result
-  double C_data[] = { -1.0500, 1.3750, 0.0833, 0.6250,
-                      0.3000, -0.7500, 0.8333, 0.7500,
-                      -0.2167, -0.0417, 0.8056, 1.3750 };
-  MatrixView C_prime(C_data, m, n, m);
-  Matrix C(n, m);
-  mat::transpose(C, C_prime);
-  mat::sub(B, C);
-  assert(mat::nrminf(B) < 1e-4);
-}

mex/sources/estimation/libmat/tests/test-qr.cc

-/*
- * 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), S2(m, n);
-  QRDecomposition QRD(m, n, m);
-
-  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);
-
-  S2 = S;
-  QRD.computeAndLeftMultByQ(S2, "N", Q);
-
-  std::cout << "Q =" << std::endl << Q << std::endl;
-
-  blas::gemm("T", "N", 1.0, Q, Q, 0.0, B);
-
-  std::cout << "Q'*Q =" << std::endl << B << std::endl;
-
-  for (size_t j = 0; j < n; j++)
-    mat::col_set(S2, j, j+1, m-j-1, 0);
-
-  std::cout << "R =" << std::endl << S2 << std::endl;
-
-  blas::gemm("N", "N", 1.0, Q, S2, 0.0, A);
-
-  std::cout << "Q*R =" << std::endl << A << std::endl;
-
-  // Check values
-  Matrix B2(m);
-  mat::set_identity(B2);
-  mat::sub(B2, B);
-  assert(mat::nrminf(B2) < 1e-4);
-
-  mat::sub(A, S);
-  assert(mat::nrminf(A) < 1e-4);
-}

mex/sources/estimation/libmat/tests/test-repmat.cc

-/*
- * 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 "Matrix.hh"
-
-int
-main(int argc, char **argv)
-{
-  size_t m = 4, n = 3;
-  Matrix S(m, n), A(5*m, 3*n), A2(5*m, 3*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 tiled:" << std::endl << S << std::endl;
-
-  Matrix   S1(m, 2 *n), S2(2*m, n);
-  std::vector<size_t>  toC(n), toR(m);
-  for (size_t i = 0; i < n; i++)
-    toC[i] = 2*i;
-
-  for (size_t j = 0; j < m; j++)
-    toR[j] = m-j-1;
-
-  mat::assignByVectors(S1, mat::nullVec, toC, S, mat::nullVec, mat::nullVec);
-  std::cout << "Matrix assigned and col reorder by vectors:" << std::endl << S1 << std::endl;
-
-  mat::assignByVectors(S2, toR, mat::nullVec, S, mat::nullVec, mat::nullVec);
-  std::cout << "Matrix assigned and row ereorder by vectors:" << std::endl << S2 << std::endl;
-
-  mat::assignByVectors(S1, toR, toC, S, mat::nullVec, mat::nullVec);
-  std::cout << "Matrix assigned by vectors:" << std::endl << S1 << std::endl;
-
-  mat::repmat(S, 5, 3, A);
-  std::cout << "Tiled=" << std::endl << A << std::endl;
-
-  A2 = A;
-  bool b = mat::isDiff(A, A2, 0.0000001);
-  std::cout << "No Diff=" << b << std::endl;
-  mat::add(A, 0.000001);
-  b = mat::isDiff(A, A2, 0.0000001);
-  std::cout << "Yes Diff=" << b << std::endl;
-  Matrix A4(5), A5(5);
-  A4.setAll(1.0);
-  A5.setAll(1.0);
-  b = mat::isDiffSym(A4, A5, 0.0000001);
-  std::cout << "No DiffSym=" << b << std::endl;
-  mat::sub(A4, 0.0001);
-  b = mat::isDiffSym(A4, A5, 0.0000001);
-  std::cout << "Yes DiffSym=" << b << std::endl;
-
-}

mex/sources/estimation/logMHMCMCposterior.cc

-/*
- * Copyright (C) 2010-2013 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 <string>
-#include <vector>
-#include <algorithm>
-#include <functional>
-
-#include "Vector.hh"
-#include "Matrix.hh"
-#include "LogPosteriorDensity.hh"
-#include "RandomWalkMetropolisHastings.hh"
-
-#include <dynmex.h>
-#if defined MATLAB_MEX_FILE
-# include "mat.h"
-#else   //  OCTAVE_MEX_FILE e.t.c.
-# include "matio.h"
-#endif
-
-#if defined(_WIN32) || defined(__CYGWIN32__) || defined(WINDOWS)
-# define DIRECTORY_SEPARATOR "\\"
-#else
-# define DIRECTORY_SEPARATOR "/"
-#endif
-
-class LogMHMCMCposteriorMexErrMsgTxtException
-{
-public:
-  std::string errMsg;
-  LogMHMCMCposteriorMexErrMsgTxtException(const std::string &msg) : errMsg(msg)
-  {
-  }
-  inline const char *getErrMsg() { return errMsg.c_str(); }
-};
-
-void
-fillEstParamsInfo(const mxArray *estim_params_info, EstimatedParameter::pType type,
-                  std::vector<EstimatedParameter> &estParamsInfo)
-{
-  const mxArray *bayestopt_ = mexGetVariablePtr("global", "bayestopt_");
-  const mxArray *bayestopt_ubp = mxGetField(bayestopt_, 0, "ub"); // upper bound
-  const mxArray *bayestopt_lbp = mxGetField(bayestopt_, 0, "lb"); // lower bound
-  const mxArray *bayestopt_p1p = mxGetField(bayestopt_, 0, "p1"); // prior mean
-  const mxArray *bayestopt_p2p = mxGetField(bayestopt_, 0, "p2"); // prior standard deviation
-  const mxArray *bayestopt_p3p = mxGetField(bayestopt_, 0, "p3"); // lower bound
-  const mxArray *bayestopt_p4p = mxGetField(bayestopt_, 0, "p4"); // upper bound
-  const mxArray *bayestopt_p6p = mxGetField(bayestopt_, 0, "p6"); // first hyper-parameter (\alpha for the BETA and GAMMA distributions, s for the INVERSE GAMMAs, expectation for the GAUSSIAN distribution, lower bound for the UNIFORM distribution).
-  const mxArray *bayestopt_p7p = mxGetField(bayestopt_, 0, "p7"); // second hyper-parameter (\beta for the BETA and GAMMA distributions, \nu for the INVERSE GAMMAs, standard deviation for the GAUSSIAN distribution, upper bound for the UNIFORM distribution).
-  const mxArray *bayestopt_jscalep = mxGetField(bayestopt_, 0, "jscale"); // MCMC jump scale
-
-  const size_t bayestopt_size = mxGetM(bayestopt_);
-  const VectorConstView bayestopt_ub(mxGetPr(bayestopt_ubp), bayestopt_size, 1);
-  const VectorConstView bayestopt_lb(mxGetPr(bayestopt_lbp), bayestopt_size, 1);
-  const VectorConstView bayestopt_p1(mxGetPr(bayestopt_p1p), bayestopt_size, 1); //=mxGetField(bayestopt_, 0, "p1");
-  const VectorConstView bayestopt_p2(mxGetPr(bayestopt_p2p), bayestopt_size, 1); //=mxGetField(bayestopt_, 0, "p2");
-  const VectorConstView bayestopt_p3(mxGetPr(bayestopt_p3p), bayestopt_size, 1); //=mxGetField(bayestopt_, 0, "p3");
-  const VectorConstView bayestopt_p4(mxGetPr(bayestopt_p4p), bayestopt_size, 1); //=mxGetField(bayestopt_, 0, "p4");
-  const VectorConstView bayestopt_p6(mxGetPr(bayestopt_p6p), bayestopt_size, 1); //=mxGetField(bayestopt_, 0, "p6");
-  const VectorConstView bayestopt_p7(mxGetPr(bayestopt_p7p), bayestopt_size, 1); //=mxGetField(bayestopt_, 0, "p7");
-  const VectorConstView bayestopt_jscale(mxGetPr(bayestopt_jscalep), bayestopt_size, 1); //=mxGetField(bayestopt_, 0, "jscale");
-
-  // loop processsing
-  size_t m = mxGetM(estim_params_info), n = mxGetN(estim_params_info);
-  MatrixConstView epi(mxGetPr(estim_params_info), m, n, m);
-  size_t bayestopt_count = estParamsInfo.size();
-
-  for (size_t i = 0; i < m; i++)
-    {
-      size_t col = 0;
-      size_t id1 = (size_t) epi(i, col++) - 1;
-      size_t id2 = 0;
-      if (type == EstimatedParameter::shock_Corr
-          || type == EstimatedParameter::measureErr_Corr)
-        id2 = (size_t) epi(i, col++) - 1;
-      col++; // Skip init_val #2 or #3
-      double par_low_bound =  bayestopt_lb(bayestopt_count); col++; //#3 epi(i, col++);
-      double par_up_bound =  bayestopt_ub(bayestopt_count); col++; //#4 epi(i, col++);
-      Prior::pShape shape = (Prior::pShape) epi(i, col++);
-      double mean = epi(i, col++);
-      double std = epi(i, col++);
-      double low_bound =  bayestopt_p3(bayestopt_count);
-      double up_bound =  bayestopt_p4(bayestopt_count);
-      double fhp =  bayestopt_p6(bayestopt_count); // double p3 = epi(i, col++);
-      double shp =  bayestopt_p7(bayestopt_count); // double p4 = epi(i, col++);
-
-      Prior *p = Prior::constructPrior(shape, mean, std, low_bound, up_bound, fhp, shp); //1.0,INFINITY);//p3, p4);
-
-      // Only one subsample
-      std::vector<size_t> subSampleIDs;
-      subSampleIDs.push_back(0);
-      estParamsInfo.push_back(EstimatedParameter(type, id1, id2, subSampleIDs,
-                                                 par_low_bound, par_up_bound, p));
-      bayestopt_count++;
-    }
-}
-
-int
-sampleMHMC(LogPosteriorDensity &lpd, RandomWalkMetropolisHastings &rwmh,
-           VectorView &steadyState, VectorConstView &estParams, VectorView &deepParams, const MatrixConstView &data,
-           MatrixView &Q, Matrix &H, size_t presampleStart, 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)
-{
-  enum {iMin, iMax};
-  int iret = 0; // return value
-  std::vector<size_t> OpenOldFile(nBlocks, 0);
-  size_t jloop = 0, irun, j; // counters
-  double dsum, dmax, dmin, sux = 0, jsux = 0;
-  std::string mhFName;
-  std::stringstream ssFName;
-#if defined MATLAB_MEX_FILE
-  MATFile *drawmat; // MCMC draws output file pointer
-  int matfStatus;
-#else   //  OCTAVE_MEX_FILE e.t.c.
-# if MATIO_MAJOR_VERSION > 1 || (MATIO_MAJOR_VERSION == 1 && MATIO_MINOR_VERSION >= 5)
-  size_t dims[2];
-  const matio_compression compression = MAT_COMPRESSION_NONE;
-# else
-  int dims[2];
-  const int compression = COMPRESSION_NONE;
-# endif
-  mat_t *drawmat;
-  matvar_t *matvar;
-  int matfStatus;
-#endif
-  FILE *fidlog;  // log file
-  size_t npar = estParams.getSize();
-  Matrix MinMax(npar, 2);
-
-  const mxArray *InitSizeArrayPtr = mexGetVariablePtr("caller", "InitSizeArray");
-  if (InitSizeArrayPtr == NULL)
-    {
-      mexPrintf("Metropolis-Hastings myinputs field InitSizeArrayPtr Initialisation failed!\n");
-      return (-1);
-    }
-  const VectorConstView InitSizeArrayVw(mxGetPr(InitSizeArrayPtr), nBlocks, 1);
-  Vector InitSizeArray(InitSizeArrayVw.getSize());
-  InitSizeArray = InitSizeArrayVw;
-  //const mxArray *flinePtr = mxGetField(myinputs, 0, "fline");
-  const mxArray *flinePtr = mexGetVariable("caller",  "fline");
-  if (flinePtr == NULL)
-    {
-      mexPrintf("Metropolis-Hastings myinputs field fline Initialisation failed!\n");
-      return (-1);
-    }
-  VectorView fline(mxGetPr(flinePtr), nBlocks, 1);
-
-  mxArray *NewFileArrayPtr = mexGetVariable("caller", "NewFile");
-  if (NewFileArrayPtr == NULL)
-    {
-      mexPrintf("Metropolis-Hastings myinputs fields NewFileArrayPtr Initialisation failed!\n");
-      return (-1);
-    }
-  VectorView NewFileVw(mxGetPr(NewFileArrayPtr), nBlocks, 1);
-  //Vector NewFile(NewFileVw.getSize());
-  //NewFile = NewFileVw;
-
-  const mxArray *MAX_nrunsPtr = mexGetVariablePtr("caller", "MAX_nruns");
-  const size_t MAX_nruns = (size_t) mxGetScalar(MAX_nrunsPtr);
-
-  const mxArray *blockStartParamsPtr = mexGetVariable("caller", "ix2");
-  MatrixView blockStartParamsMxVw(mxGetPr(blockStartParamsPtr), nBlocks, npar, nBlocks);
-  Vector startParams(npar);
-
-  const mxArray *mxFirstLogLikPtr = mexGetVariable("caller", "ilogpo2");
-  VectorView FirstLogLiK(mxGetPr(mxFirstLogLikPtr), nBlocks, 1);
-
-  const mxArray *record = mexGetVariable("caller", "record");
-  //const mxArray *record = mxGetField(myinputs, 0, "record");
-  if (record == NULL)
-    {
-      mexPrintf("Metropolis-Hastings record Initialisation failed!\n");
-      return (-1);
-    }
-  mxArray *AcceptationRatesPtr = mxGetField(record, 0, "AcceptationRates");
-  if (AcceptationRatesPtr == NULL)
-    {
-      mexPrintf("Metropolis-Hastings record AcceptationRatesPtr Initialisation failed!\n");
-      return (-1);
-    }
-  VectorView AcceptationRates(mxGetPr(AcceptationRatesPtr), nBlocks, 1);
-
-  mxArray *mxLastParametersPtr = mxGetField(record, 0, "LastParameters");
-  MatrixView LastParameters(mxGetPr(mxLastParametersPtr), nBlocks, npar, nBlocks);
-  LastParameters = blockStartParamsMxVw;
-
-  mxArray *mxLastLogLikPtr = mxGetField(record, 0, "LastLogLiK");
-  VectorView LastLogLiK(mxGetPr(mxLastLogLikPtr), nBlocks, 1);
-
-  mxArray *mxMhLogPostDensPtr = 0;
-  mxArray *mxMhParamDrawsPtr = 0;
-  size_t currInitSizeArray = 0;
-
-#if defined MATLAB_MEX_FILE
-  // Waitbar
-  mxArray *waitBarRhs[3], *waitBarLhs[1];
-  waitBarRhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
-  std::string barTitle;
-  std::stringstream ssbarTitle;
-  if (console_mode == 0)
-    {
-      ssbarTitle.clear();
-      ssbarTitle.str("");
-      ssbarTitle << "Please wait... Metropolis-Hastings  " << fblock << "/" << nBlocks << "  ...";
-      barTitle = ssbarTitle.str();
-      waitBarRhs[1] = mxCreateString(barTitle.c_str());
-      *mxGetPr(waitBarRhs[0]) = (double) 0.0;
-      mexCallMATLAB(1, waitBarLhs, 2, waitBarRhs, "waitbar");
-      if (waitBarRhs[1])
-        mxDestroyArray(waitBarRhs[1]);
-      waitBarRhs[1] = waitBarLhs[0];
-    }
-#endif
-
-  for (size_t b = fblock; b <= nBlocks; ++b)
-    {
-      jloop = jloop+1;
-
-#if defined MATLAB_MEX_FILE
-      if ((load_mh_file != 0)  && (fline(b) > 1) && OpenOldFile[b])
-        {
-          //  load(['./' MhDirectoryName '/' ModelName '_mh' int2str(NewFile(b)) '_blck' int2str(b) '.mat'])
-          ssFName.clear();
-          ssFName.str("");
-          ssFName << resultsFileStem << DIRECTORY_SEPARATOR << "metropolis" << DIRECTORY_SEPARATOR << resultsFileStem << "_mh" << (size_t) NewFileVw(b-1) << "_blck" << b << ".mat";
-          mhFName = ssFName.str();
-          drawmat = matOpen(mhFName.c_str(), "r");
-          mexPrintf("MHMCMC: Using interim partial draws file %s \n", mhFName.c_str());
-          if (drawmat == 0)
-            {
-              fline(b) = 1;
-              mexPrintf("Error in MH: Can not open old draws Mat file for reading:  %s \n  \
-                  Starting a new file instead! \n", mhFName.c_str());
-            }
-          else
-            {
-              currInitSizeArray = (size_t) InitSizeArray(b-1);
-              mxMhParamDrawsPtr = matGetVariable(drawmat, "x2");
-              mxMhLogPostDensPtr = matGetVariable(drawmat, "logpo2");
-              matClose(drawmat);
-              OpenOldFile[b] = 1;
-            }
-        } // end if
-      if (console_mode == 0)
-        {
-          ssbarTitle.clear();
-          ssbarTitle.str("");
-          ssbarTitle << "Please wait... Metropolis-Hastings  " << b << "/" << nBlocks << "  ...";
-          barTitle = ssbarTitle.str();
-          waitBarRhs[2] = mxCreateString(barTitle.c_str());
-          //strcpy( *mxGetPr(waitBarRhs[1]), mhFName.c_str());
-          *mxGetPr(waitBarRhs[0]) = (double) 0.0;
-          mexCallMATLAB(0, NULL, 3, waitBarRhs, "waitbar");
-          mxDestroyArray(waitBarRhs[2]);
-        }
-#else //if defined OCTAVE_MEX_FILE
-      if ((load_mh_file != 0)  && (fline(b) > 1) && OpenOldFile[b])
-        {
-          //  load(['./' MhDirectoryName '/' ModelName '_mh' int2str(NewFile(b)) '_blck' int2str(b) '.mat'])
-          if ((currInitSizeArray != (size_t) InitSizeArray(b-1)) &&   OpenOldFile[b] != 1)
-            {
-              // new or different size result arrays/matrices
-              currInitSizeArray = (size_t) InitSizeArray(b-1);
-              if (mxMhLogPostDensPtr)
-                mxDestroyArray(mxMhLogPostDensPtr);                                                                                                          // log post density array
-              mxMhLogPostDensPtr = mxCreateDoubleMatrix(currInitSizeArray, 1, mxREAL);
-              if (mxMhLogPostDensPtr == NULL)
-                {
-                  mexPrintf("Metropolis-Hastings mxMhLogPostDensPtr Initialisation failed!\n");
-                  return (-1);
-                }
-              if (mxMhParamDrawsPtr)
-                mxDestroyArray(mxMhParamDrawsPtr);                                                                                                        // accepted MCMC MH draws
-              mxMhParamDrawsPtr =  mxCreateDoubleMatrix(currInitSizeArray, npar,  mxREAL);
-              if (mxMhParamDrawsPtr == NULL)
-                {
-                  mexPrintf("Metropolis-Hastings mxMhParamDrawsPtr Initialisation failed!\n");
-                  return (-1);
-                }
-            }
-
-          ssFName.clear();
-          ssFName.str("");
-          ssFName << resultsFileStem << DIRECTORY_SEPARATOR << "metropolis" << DIRECTORY_SEPARATOR << resultsFileStem << "_mh" << (size_t) NewFileVw(b-1) << "_blck" << b << ".mat";
-          mhFName = ssFName.str();
-          drawmat = Mat_Open(mhFName.c_str(), MAT_ACC_RDONLY);
-          if (drawmat == NULL)
-            {
-              fline(b) = 1;
-              mexPrintf("Error in MH: Can not open old draws Mat file for reading:  %s \n  \
-                  Starting a new file instead! \n", mhFName.c_str());
-            }
-          else
-            {
-              int start[2] = {0, 0}, edge[2], stride[2] = {1, 1}, err = 0;
-              mexPrintf("MHMCMC: Using interim partial draws file %s \n", mhFName.c_str());
-              //              matvar = Mat_VarReadInfo(drawmat, "x2");
-              matvar = Mat_VarReadInfo(drawmat, (char *) "x2");
-              if (matvar == NULL)
-                {
-                  fline(b) = 1;
-                  mexPrintf("Error in MH: Can not read old draws Mat file for reading:  %s \n  \
-                      Starting a new file instead! \n", mhFName.c_str());
-                }
-              else
-                {
-                  // GetVariable(drawmat, "x2");
-                  edge[0] = matvar->dims[0];
-                  edge[1] = matvar->dims[1];
-                  err = Mat_VarReadData(drawmat, matvar, mxGetPr(mxMhParamDrawsPtr), start, stride, edge);
-                  if (err)
-                    {
-                      fline(b) = 1;
-                      mexPrintf("Error in MH: Can not retreive old draws from Mat file:  %s \n  \
-                          Starting a new file instead! \n", mhFName.c_str());
-                    }
-                  Mat_VarFree(matvar);
-                }
-              //mxMhLogPostDensPtr = Mat_GetVariable(drawmat, "logpo2");
-              matvar = Mat_VarReadInfo(drawmat, (char *) "logpo2");
-              if (matvar == NULL)
-                {
-                  fline(b) = 1;
-                  mexPrintf("Error in MH: Can not read old logPos Mat file for reading:  %s \n  \
-                      Starting a new file instead! \n", mhFName.c_str());
-                }
-              else
-                {
-                  // GetVariable(drawmat, "x2");
-                  edge[0] = matvar->dims[0];
-                  edge[1] = matvar->dims[1];
-                  err = Mat_VarReadData(drawmat, matvar, mxGetPr(mxMhLogPostDensPtr), start, stride, edge);
-                  if (err)
-                    {
-                      fline(b) = 1;
-                      mexPrintf("Error in MH: Can not retreive old logPos from Mat file:  %s \n  \
-                          Starting a new file instead! \n", mhFName.c_str());
-                    }
-                  Mat_VarFree(matvar);
-                }
-              Mat_Close(drawmat);
-              OpenOldFile[b] = 1;
-            }
-        } // end if
-
-#endif
-
-      VectorView LastParametersRow = mat::get_row(LastParameters, b-1);
-
-      sux = 0.0;
-      jsux = 0;
-      irun = (size_t) fline(b-1);
-      j = 0; //1;
-      while (j < nruns(b-1))
-        {
-          if ((currInitSizeArray != (size_t) InitSizeArray(b-1)) &&   OpenOldFile[b] != 1)
-            {
-              // new or different size result arrays/matrices
-              currInitSizeArray = (size_t) InitSizeArray(b-1);
-              if (mxMhLogPostDensPtr)
-                mxDestroyArray(mxMhLogPostDensPtr);                                                                                                          // log post density array
-              mxMhLogPostDensPtr = mxCreateDoubleMatrix(currInitSizeArray, 1, mxREAL);
-              if (mxMhLogPostDensPtr == NULL)
-                {
-                  mexPrintf("Metropolis-Hastings mxMhLogPostDensPtr Initialisation failed!\n");
-                  return (-1);
-                }
-              if (mxMhParamDrawsPtr)
-                mxDestroyArray(mxMhParamDrawsPtr);                                                                                                        // accepted MCMC MH draws
-              mxMhParamDrawsPtr =  mxCreateDoubleMatrix(currInitSizeArray, npar,  mxREAL);
-              if (mxMhParamDrawsPtr == NULL)
-                {
-                  mexPrintf("Metropolis-Hastings mxMhParamDrawsPtr Initialisation failed!\n");
-                  return (-1);
-                }
-            }
-          startParams = LastParametersRow;
-          VectorView mhLogPostDens(mxGetPr(mxMhLogPostDensPtr), currInitSizeArray, (size_t) 1);
-          MatrixView mhParamDraws(mxGetPr(mxMhParamDrawsPtr), currInitSizeArray, npar, currInitSizeArray);
-          try
-            {
-              jsux = rwmh.compute(mhLogPostDens, mhParamDraws, steadyState, startParams, deepParams, data, Q, H,
-                                  presampleStart, irun, currInitSizeArray, lpd, pdd, epd);
-              irun = currInitSizeArray;
-              sux += jsux*currInitSizeArray;
-              j += currInitSizeArray; //j=j+1;
-            }
-          catch (const TSException &tse)
-            {
-              iret = -100;
-              mexPrintf(" TSException Exception in RandomWalkMH dynamic_dll: %s \n", (tse.getMessage()).c_str());
-              goto cleanup;
-            }
-          catch (const DecisionRules::BlanchardKahnException &bke)
-            {
-              iret = -90;
-              mexPrintf(" Too many Blanchard-Kahn Exceptions in RandomWalkMH : n_fwrd_vars %d n_explosive_eigenvals %d \n", bke.n_fwrd_vars, bke.n_explosive_eigenvals);
-              goto cleanup;
-            }
-          catch (const GeneralizedSchurDecomposition::GSDException &gsde)
-            {
-              iret = -80;
-              mexPrintf(" GeneralizedSchurDecomposition Exception in RandomWalkMH: info %d, n %d  \n", gsde.info, gsde.n);
-              goto cleanup;
-            }
-          catch (const LUSolver::LUException &lue)
-            {
-              iret = -70;
-              mexPrintf(" LU Exception in RandomWalkMH : info %d \n", lue.info);
-              goto cleanup;
-            }
-          catch (const VDVEigDecomposition::VDVEigException &vdve)
-            {
-              iret = -60;
-              mexPrintf(" VDV Eig Exception in RandomWalkMH : %s ,  info: %d\n", vdve.message.c_str(), vdve.info);
-              goto cleanup;
-            }
-          catch (const DiscLyapFast::DLPException &dlpe)
-            {
-              iret = -50;
-              mexPrintf(" Lyapunov solver Exception in RandomWalkMH : %s ,  info: %d\n", dlpe.message.c_str(), dlpe.info);
-              goto cleanup;
-            }
-          catch (const std::runtime_error &re)
-            {
-              iret = -3;
-              mexPrintf(" Runtime Error Exception in RandomWalkMH: %s \n", re.what());
-              goto cleanup;
-            }
-          catch (const std::exception &e)
-            {
-              iret = -2;
-              mexPrintf(" Standard System Exception in RandomWalkMH: %s \n", e.what());
-              goto cleanup;
-            }
-          catch (...)
-            {
-              iret = -1000;
-              mexPrintf(" Unknown unhandled Exception in RandomWalkMH! %s \n");
-              goto cleanup;
-            }
-
-#if defined MATLAB_MEX_FILE
-          if (console_mode)
-            mexPrintf("   MH: Computing Metropolis-Hastings (chain %d/%d): %3.f \b%% done, acceptance rate: %3.f \b%%\r", b, nBlocks, 100 * j/nruns(b-1), 100 * sux / j);
-          else
-            {
-              // Waitbar
-              ssbarTitle.clear();
-              ssbarTitle.str("");
-              ssbarTitle << "Metropolis-Hastings : " << b << "/" << nBlocks << " Acceptance: " << 100 * sux/j << "%";
-              barTitle = ssbarTitle.str();
-              waitBarRhs[2] = mxCreateString(barTitle.c_str());
-              *mxGetPr(waitBarRhs[0]) = j / nruns(b-1);
-              mexCallMATLAB(0, NULL, 3, waitBarRhs, "waitbar");
-              mxDestroyArray(waitBarRhs[2]);
-
-            }
-
-          // % Now I save the simulations
-          // save draw  2 mat file ([MhDirectoryName '/' ModelName '_mh' int2str(NewFile(b)) '_blck' int2str(b) '.mat'],'x2','logpo2');
-          ssFName.clear();
-          ssFName.str("");
-          ssFName << resultsFileStem << DIRECTORY_SEPARATOR << "metropolis" << DIRECTORY_SEPARATOR << resultsFileStem << "_mh" << (size_t) NewFileVw(b-1) << "_blck" << b << ".mat";
-          mhFName = ssFName.str();
-          drawmat = matOpen(mhFName.c_str(), "w");
-          if (drawmat == 0)
-            {
-              mexPrintf("Error in MH: Can not open draws Mat file for writing:  %s \n", mhFName.c_str());
-              exit(1);
-            }
-          matfStatus = matPutVariable(drawmat, "x2", mxMhParamDrawsPtr);
-          if (matfStatus)
-            {
-              mexPrintf("Error in MH: Can not use draws Mat file for writing:  %s \n", mhFName.c_str());
-              exit(1);
-            }
-          matfStatus = matPutVariable(drawmat, "logpo2", mxMhLogPostDensPtr);
-          if (matfStatus)
-            {
-              mexPrintf("Error in MH: Can not usee draws Mat file for writing:  %s \n", mhFName.c_str());
-              exit(1);
-            }
-          matClose(drawmat);
-#else
-
-          printf("   MH: Computing Metropolis-Hastings (chain %ld/%ld): %3.f \b%% done, acceptance rate: %3.f \b%%\r", b, nBlocks, 100 * j/nruns(b-1), 100 * sux / j);
-          // % Now I save the simulations
-          // save draw  2 mat file ([MhDirectoryName '/' ModelName '_mh' int2str(NewFile(b)) '_blck' int2str(b) '.mat'],'x2','logpo2');
-          ssFName.clear();
-          ssFName.str("");
-          ssFName << resultsFileStem << DIRECTORY_SEPARATOR << "metropolis" << DIRECTORY_SEPARATOR << resultsFileStem << "_mh" << (size_t) NewFileVw(b-1) << "_blck" << b << ".mat";
-          mhFName = ssFName.str();
-
-          drawmat = Mat_Open(mhFName.c_str(), MAT_ACC_RDWR);
-          if (drawmat == 0)
-            {
-              mexPrintf("Error in MH: Can not open draws Mat file for writing:  %s \n", mhFName.c_str());
-              exit(1);
-            }
-          dims[0] = currInitSizeArray;
-          dims[1] = npar;
-          matvar = Mat_VarCreate("x2", MAT_C_DOUBLE, MAT_T_DOUBLE, 2, dims, mxGetPr(mxMhParamDrawsPtr), 0);
-          matfStatus = Mat_VarWrite(drawmat, matvar, compression);
-          Mat_VarFree(matvar);
-          if (matfStatus)
-            {
-              mexPrintf("Error in MH: Can not use draws Mat file for writing:  %s \n", mhFName.c_str());
-              exit(1);
-            }
-          //matfStatus = matPutVariable(drawmat, "logpo2", mxMhLogPostDensPtr);
-          dims[1] = 1;
-          matvar = Mat_VarCreate("logpo2", MAT_C_DOUBLE, MAT_T_DOUBLE, 2, dims, mxGetPr(mxMhLogPostDensPtr), 0);
-          matfStatus = Mat_VarWrite(drawmat, matvar, compression);
-          Mat_VarFree(matvar);
-          if (matfStatus)
-            {
-              mexPrintf("Error in MH: Can not usee draws Mat file for writing:  %s \n", mhFName.c_str());
-              exit(1);
-            }
-          Mat_Close(drawmat);
-#endif
-
-          // save log to fidlog = fopen([MhDirectoryName '/metropolis.log'],'a');
-          ssFName.str("");
-          ssFName << resultsFileStem << DIRECTORY_SEPARATOR << "metropolis" << DIRECTORY_SEPARATOR << "metropolis.log";
-          mhFName = ssFName.str();
-          fidlog = fopen(mhFName.c_str(), "a");
-          fprintf(fidlog, "\n");
-          fprintf(fidlog, "%% Mh%dBlck%lu ( %s %s )\n", (int) NewFileVw(b-1), b, __DATE__, __TIME__);
-          fprintf(fidlog, " \n");
-          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: %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: %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: %lu : %f \n", i+1, MinMax(i, iMax));
-          fprintf(fidlog, "    log2po: %f \n", dmax);
-          fprintf(fidlog, " \n");
-          fclose(fidlog);
-
-          jsux = 0;
-          LastParametersRow = mat::get_row(mhParamDraws, currInitSizeArray-1); //x2(end,:);
-          LastLogLiK(b-1) = mhLogPostDens(currInitSizeArray-1); //logpo2(end);
-          InitSizeArray(b-1) = std::min((size_t) nruns(b-1)-j, MAX_nruns);
-          // initialization of next file if necessary
-          if (InitSizeArray(b-1))
-            {
-              NewFileVw(b-1)++; // = NewFile(b-1) + 1;
-              irun = 1;
-            } // end
-          //irun++;
-        } // end  while % End of the simulations for one mh-block.
-      //record.
-      AcceptationRates(b-1) = sux/j;
-      OpenOldFile[b] = 0;
-    } // end % End of the loop over the mh-blocks.
-
-  if (mexPutVariable("caller", "record_AcceptationRates", AcceptationRatesPtr))
-    mexPrintf("MH Warning: due to error record_AcceptationRates is NOT set !! \n");
-
-  if (mexPutVariable("caller", "record_LastParameters", mxLastParametersPtr))
-    mexPrintf("MH Warning: due to error record_MhParamDraw is NOT set !! \n");
-
-  if (mexPutVariable("caller", "record_LastLogLiK", mxLastLogLikPtr))
-    mexPrintf("MH Warning: due to error record_LastLogLiK is NOT set !! \n");
-
-  //NewFileVw = NewFile;
-  if (mexPutVariable("caller", "NewFile", NewFileArrayPtr))
-    mexPrintf("MH Warning: due to error NewFile is NOT set !! \n");
-
-  // Cleanup
-  mexPrintf("MH Cleanup !! \n");
-
- cleanup:
-  if (mxMhLogPostDensPtr)
-    mxDestroyArray(mxMhLogPostDensPtr);                                            // delete log post density array
-  if (mxMhParamDrawsPtr)
-    mxDestroyArray(mxMhParamDrawsPtr);                                            // delete accepted MCMC MH draws
-
-#ifdef MATLAB_MEX_FILE
-  // Waitbar
-  if (console_mode == 0)
-    {
-      // Bellow call to close waitbar seems to cause crashes and it is for
-      // now left commented out and the waitbar neeeds to be closed manually
-      // alternativelly, call with options_.console_mode=1;
-      //mexCallMATLAB(0, NULL, 1, waitBarLhs, "close");
-      //mxDestroyArray(waitBarLhs[0]);
-      mxDestroyArray(waitBarRhs[1]);
-      mxDestroyArray(waitBarRhs[0]);
-    }
-#endif
-
-  // return error code or last line run in the last MH block sub-array
-  if (iret == 0)
-    iret = (int) irun;
-  return iret;
-
-}
-
-int
-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 *options_ = mexGetVariablePtr("global", "options_");
-  const mxArray *estim_params_ = mexGetVariablePtr("global", "estim_params_");
-  const mxArray *bayestopt_ = mexGetVariablePtr("global", "bayestopt_");
-
-  double loglinear = *mxGetPr(mxGetField(options_, 0, "loglinear"));
-  if (loglinear == 1)
-    throw LogMHMCMCposteriorMexErrMsgTxtException("Option loglinear is not supported");
-
-  if (*mxGetPr(mxGetField(options_, 0, "endogenous_prior")) == 1)
-    throw LogMHMCMCposteriorMexErrMsgTxtException("Option endogenous_prior is not supported");
-
-  double with_trend = *mxGetPr(mxGetField(bayestopt_, 0, "with_trend"));
-  if (with_trend == 1)
-    throw LogMHMCMCposteriorMexErrMsgTxtException("Observation trends are not supported");
-
-  // Construct arguments of constructor of LogLikelihoodMain
-  char *fName = mxArrayToString(mxGetField(M_, 0, "fname"));
-  std::string resultsFileStem(fName);
-  std::string basename(fName);
-  mxFree(fName);
-
-  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_estParams = estParams.getSize();
-
-  std::vector<size_t> zeta_fwrd, zeta_back, zeta_mixed, zeta_static;
-  const mxArray *lli_mx = mxGetField(M_, 0, "lead_lag_incidence");
-  MatrixConstView lli(mxGetPr(lli_mx), mxGetM(lli_mx), mxGetN(lli_mx), mxGetM(lli_mx));
-  if (lli.getRows() != 3 || lli.getCols() != n_endo)
-    throw LogMHMCMCposteriorMexErrMsgTxtException("Error in logMCMCposterior: Incorrect lead/lag incidence matrix");
-  for (size_t i = 0; i < n_endo; i++)
-    {
-      if (lli(0, i) == 0 && lli(2, i) == 0)
-        zeta_static.push_back(i);
-      else if (lli(0, i) != 0 && lli(2, i) == 0)
-        zeta_back.push_back(i);
-      else if (lli(0, i) == 0 && lli(2, i) != 0)
-        zeta_fwrd.push_back(i);
-      else
-        zeta_mixed.push_back(i);
-    }
-
-  double qz_criterium = *mxGetPr(mxGetField(options_, 0, "qz_criterium"));
-  double lyapunov_tol = *mxGetPr(mxGetField(options_, 0, "lyapunov_complex_threshold"));
-  double riccati_tol = *mxGetPr(mxGetField(options_, 0, "riccati_tol"));
-  size_t presample = (size_t) *mxGetPr(mxGetField(options_, 0, "presample"));
-  size_t console_mode = (size_t) *mxGetPr(mxGetField(options_, 0, "console_mode"));
-  size_t load_mh_file = (size_t) *mxGetPr(mxGetField(options_, 0, "load_mh_file"));
-
-  std::vector<size_t> varobs;
-  const mxArray *varobs_mx = mxGetField(options_, 0, "varobs_id");
-  if (mxGetM(varobs_mx) != 1)
-    throw LogMHMCMCposteriorMexErrMsgTxtException("Error in logMCMCposterior: options_.varobs_id must be a row vector");
-  size_t n_varobs = mxGetN(varobs_mx);
-  std::transform(mxGetPr(varobs_mx), mxGetPr(varobs_mx) + n_varobs, back_inserter(varobs),
-                 std::bind2nd(std::minus<size_t>(), 1));
-
-  if (data.getRows() != n_varobs)
-    throw LogMHMCMCposteriorMexErrMsgTxtException("Error in logMCMCposterior: Data does not have as many rows as there are observed variables");
-
-  std::vector<EstimationSubsample> estSubsamples;
-  estSubsamples.push_back(EstimationSubsample(0, data.getCols() - 1));
-
-  std::vector<EstimatedParameter> estParamsInfo;
-  fillEstParamsInfo(mxGetField(estim_params_, 0, "var_exo"), EstimatedParameter::shock_SD,
-                    estParamsInfo);
-  fillEstParamsInfo(mxGetField(estim_params_, 0, "var_endo"), EstimatedParameter::measureErr_SD,
-                    estParamsInfo);
-  fillEstParamsInfo(mxGetField(estim_params_, 0, "corrx"), EstimatedParameter::shock_Corr,
-                    estParamsInfo);
-  fillEstParamsInfo(mxGetField(estim_params_, 0, "corrn"), EstimatedParameter::measureErr_Corr,
-                    estParamsInfo);
-  fillEstParamsInfo(mxGetField(estim_params_, 0, "param_vals"), EstimatedParameter::deepPar,
-                    estParamsInfo);
-  EstimatedParametersDescription epd(estSubsamples, estParamsInfo);
-
-
-  bool noconstant = (bool) *mxGetPr(mxGetField(options_, 0, "noconstant"));
-
-  // Allocate LogPosteriorDensity object
-  LogPosteriorDensity lpd(basename, epd, n_endo, n_exo, zeta_fwrd, zeta_back, zeta_mixed, zeta_static,
-                          qz_criterium, varobs, riccati_tol, lyapunov_tol, noconstant);
-
-  // Construct MHMCMC Sampler
-  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);
-  const Matrix Jscale(n_estParams);
-  const VectorConstView vJscale(mxGetPr(mxGetField(bayestopt_, 0, "jscale")), n_estParams, 1);
-  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, estParams, deepParams, data, Q, H, presample,
-                                        nMHruns, fblock, nBlocks, pdd, epd, resultsFileStem, console_mode, load_mh_file);
-
-  // Cleanups
-  for (std::vector<EstimatedParameter>::iterator it = estParamsInfo.begin();
-       it != estParamsInfo.end(); it++)
-    delete it->prior;
-
-  return lastMHblockArrayLine;
-}
-
-void
-mexFunction(int nlhs, mxArray *plhs[],
-            int nrhs, const mxArray *prhs[])
-{
-  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.");
-
-  // Check and retrieve the 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");
-
-  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);
-
-  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, steadyState, deepParams, Q, H);
-      plhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL);
-      *mxGetPr(plhs[1]) = (double) lastMHblockArrayLine;
-    }
-  catch (LogMHMCMCposteriorMexErrMsgTxtException e)
-    {
-      DYN_MEX_FUNC_ERR_MSG_TXT(e.getErrMsg());
-    }
-  plhs[0] = mxCreateDoubleScalar(0);
-}

mex/sources/estimation/logposterior.cc

-/*
- * Copyright (C) 2010-2013 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 <string>
-#include <vector>
-#include <algorithm>
-#include <functional>
-#include <sstream>
-
-#include "Vector.hh"
-#include "Matrix.hh"
-#include "LogPosteriorDensity.hh"
-
-#include <dynmex.h>
-
-class LogposteriorMexErrMsgTxtException
-{
-public:
-  std::string errMsg;
-  LogposteriorMexErrMsgTxtException(const std::string &msg) : errMsg(msg)
-  {
-  }
-  inline const char *getErrMsg() { return errMsg.c_str(); }
-};
-
-void
-fillEstParamsInfo(const mxArray *bayestopt_, const mxArray *estim_params_info, EstimatedParameter::pType type,
-                  std::vector<EstimatedParameter> &estParamsInfo)
-{
-  const mxArray *bayestopt_ubp = mxGetField(bayestopt_, 0, "ub"); // upper bound
-  const mxArray *bayestopt_lbp = mxGetField(bayestopt_, 0, "lb"); // lower bound
-  const mxArray *bayestopt_p1p = mxGetField(bayestopt_, 0, "p1"); // prior mean
-  const mxArray *bayestopt_p2p = mxGetField(bayestopt_, 0, "p2"); // prior standard deviation
-  const mxArray *bayestopt_p3p = mxGetField(bayestopt_, 0, "p3"); // lower bound
-  const mxArray *bayestopt_p4p = mxGetField(bayestopt_, 0, "p4"); // upper bound
-  const mxArray *bayestopt_p6p = mxGetField(bayestopt_, 0, "p6"); // first hyper-parameter (\alpha for the BETA and GAMMA distributions, s for the INVERSE GAMMAs, expectation for the GAUSSIAN distribution, lower bound for the UNIFORM distribution).
-  const mxArray *bayestopt_p7p = mxGetField(bayestopt_, 0, "p7"); // second hyper-parameter (\beta for the BETA and GAMMA distributions, \nu for the INVERSE GAMMAs, standard deviation for the GAUSSIAN distribution, upper bound for the UNIFORM distribution).
-  const mxArray *bayestopt_jscalep = mxGetField(bayestopt_, 0, "jscale"); // MCMC jump scale
-
-  const size_t bayestopt_size = mxGetM(bayestopt_);
-  const VectorConstView bayestopt_ub(mxGetPr(bayestopt_ubp), bayestopt_size, 1);
-  const VectorConstView bayestopt_lb(mxGetPr(bayestopt_lbp), bayestopt_size, 1);
-  const VectorConstView bayestopt_p1(mxGetPr(bayestopt_p1p), bayestopt_size, 1); //=mxGetField(bayestopt_, 0, "p1");
-  const VectorConstView bayestopt_p2(mxGetPr(bayestopt_p2p), bayestopt_size, 1); //=mxGetField(bayestopt_, 0, "p2");
-  const VectorConstView bayestopt_p3(mxGetPr(bayestopt_p3p), bayestopt_size, 1); //=mxGetField(bayestopt_, 0, "p3");
-  const VectorConstView bayestopt_p4(mxGetPr(bayestopt_p4p), bayestopt_size, 1); //=mxGetField(bayestopt_, 0, "p4");
-  const VectorConstView bayestopt_p6(mxGetPr(bayestopt_p6p), bayestopt_size, 1); //=mxGetField(bayestopt_, 0, "p6");
-  const VectorConstView bayestopt_p7(mxGetPr(bayestopt_p7p), bayestopt_size, 1); //=mxGetField(bayestopt_, 0, "p7");
-  const VectorConstView bayestopt_jscale(mxGetPr(bayestopt_jscalep), bayestopt_size, 1); //=mxGetField(bayestopt_, 0, "jscale");
-
-  // loop processsing
-  size_t m = mxGetM(estim_params_info), n = mxGetN(estim_params_info);
-  MatrixConstView epi(mxGetPr(estim_params_info), m, n, m);
-  size_t bayestopt_count = estParamsInfo.size();
-
-  for (size_t i = 0; i < m; i++)
-    {
-      size_t col = 0;
-      size_t id1 = (size_t) epi(i, col++) - 1;
-      size_t id2 = 0;
-      if (type == EstimatedParameter::shock_Corr
-          || type == EstimatedParameter::measureErr_Corr)
-        id2 = (size_t) epi(i, col++) - 1;
-      col++; // Skip init_val #2 or #3
-      double par_low_bound =  bayestopt_lb(bayestopt_count); col++; //#3 epi(i, col++);
-      double par_up_bound =  bayestopt_ub(bayestopt_count); col++; //#4 epi(i, col++);
-      Prior::pShape shape = (Prior::pShape) epi(i, col++);
-      double mean = epi(i, col++);
-      double std = epi(i, col++);
-      double low_bound =  bayestopt_p3(bayestopt_count);
-      double up_bound =  bayestopt_p4(bayestopt_count);
-      double fhp =  bayestopt_p6(bayestopt_count); // double p3 = epi(i, col++);
-      double shp =  bayestopt_p7(bayestopt_count); // double p4 = epi(i, col++);
-
-      Prior *p = Prior::constructPrior(shape, mean, std, low_bound, up_bound, fhp, shp); //1.0,INFINITY);//p3, p4);
-
-      // Only one subsample
-      std::vector<size_t> subSampleIDs;
-      subSampleIDs.push_back(0);
-      estParamsInfo.push_back(EstimatedParameter(type, id1, id2, subSampleIDs,
-                                                 par_low_bound, par_up_bound, p));
-      bayestopt_count++;
-    }
-}
-
-template <class VEC1, class VEC2>
-double
-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,
-	     VectorView &deepParams, Matrix &H, MatrixView &Q)
-{
-  double loglinear = *mxGetPr(mxGetField(options_, 0, "loglinear"));
-  if (loglinear == 1)
-    throw LogposteriorMexErrMsgTxtException("Option loglinear is not supported");
-
-  if (*mxGetPr(mxGetField(options_, 0, "endogenous_prior")) == 1)
-    throw LogposteriorMexErrMsgTxtException("Option endogenous_prior is not supported");
-
-  double with_trend = *mxGetPr(mxGetField(bayestopt_, 0, "with_trend"));
-  if (with_trend == 1)
-    throw LogposteriorMexErrMsgTxtException("Observation trends are not supported");
-
-  // Construct arguments of constructor of LogLikelihoodMain
-  char *fName = mxArrayToString(mxGetField(M_, 0, "fname"));
-  std::string basename(fName);
-  mxFree(fName);
-
-  size_t n_endo = (size_t) *mxGetPr(mxGetField(M_, 0, "endo_nbr"));
-  size_t n_exo = (size_t) *mxGetPr(mxGetField(M_, 0, "exo_nbr"));
-
-  std::vector<size_t> zeta_fwrd, zeta_back, zeta_mixed, zeta_static;
-  const mxArray *lli_mx = mxGetField(M_, 0, "lead_lag_incidence");
-  MatrixConstView lli(mxGetPr(lli_mx), mxGetM(lli_mx), mxGetN(lli_mx), mxGetM(lli_mx));
-  if (lli.getRows() != 3)
-    throw LogposteriorMexErrMsgTxtException("Purely backward or purely forward models are not supported");
-  if (lli.getCols() != n_endo)
-    throw LogposteriorMexErrMsgTxtException("Incorrect lead/lag incidence matrix");
-
-  for (size_t i = 0; i < n_endo; i++)
-    {
-      if (lli(0, i) == 0 && lli(2, i) == 0)
-        zeta_static.push_back(i);
-      else if (lli(0, i) != 0 && lli(2, i) == 0)
-        zeta_back.push_back(i);
-      else if (lli(0, i) == 0 && lli(2, i) != 0)
-        zeta_fwrd.push_back(i);
-      else
-        zeta_mixed.push_back(i);
-    }
-
-  double qz_criterium = *mxGetPr(mxGetField(options_, 0, "qz_criterium"));
-  double lyapunov_tol = *mxGetPr(mxGetField(options_, 0, "lyapunov_complex_threshold"));
-  double riccati_tol = *mxGetPr(mxGetField(options_, 0, "riccati_tol"));
-  size_t presample = (size_t) *mxGetPr(mxGetField(options_, 0, "presample"));
-
-  std::vector<size_t> varobs;
-  const mxArray *varobs_mx = mxGetField(options_, 0, "varobs_id");
-  if (mxGetM(varobs_mx) != 1)
-    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));
-
-  if (data.getRows() != n_varobs)
-    throw LogposteriorMexErrMsgTxtException("Data does not have as many rows as there are observed variables");
-
-  std::vector<EstimationSubsample> estSubsamples;
-  estSubsamples.push_back(EstimationSubsample(0, data.getCols() - 1));
-
-  std::vector<EstimatedParameter> estParamsInfo;
-  fillEstParamsInfo(bayestopt_, mxGetField(estim_params_, 0, "var_exo"), EstimatedParameter::shock_SD,
-                    estParamsInfo);
-  fillEstParamsInfo(bayestopt_, mxGetField(estim_params_, 0, "var_endo"), EstimatedParameter::measureErr_SD,
-                    estParamsInfo);
-  fillEstParamsInfo(bayestopt_, mxGetField(estim_params_, 0, "corrx"), EstimatedParameter::shock_Corr,
-                    estParamsInfo);
-  fillEstParamsInfo(bayestopt_, mxGetField(estim_params_, 0, "corrn"), EstimatedParameter::measureErr_Corr,
-                    estParamsInfo);
-  fillEstParamsInfo(bayestopt_, mxGetField(estim_params_, 0, "param_vals"), EstimatedParameter::deepPar,
-                    estParamsInfo);
-
-  EstimatedParametersDescription epd(estSubsamples, estParamsInfo);
-
-  bool noconstant = (bool) *mxGetPr(mxGetField(options_, 0, "noconstant"));
-
-  // Allocate LogPosteriorDensity object
-  LogPosteriorDensity lpd(basename, epd, n_endo, n_exo, zeta_fwrd, zeta_back, zeta_mixed, zeta_static,
-                          qz_criterium, varobs, riccati_tol, lyapunov_tol, noconstant);
-
-  // Construct arguments of compute() method
-
-  // Compute the posterior
-  double logPD = lpd.compute(steadyState, estParams, deepParams, data, Q, H, presample);
-
-  // Cleanups
-  for (std::vector<EstimatedParameter>::iterator it = estParamsInfo.begin();
-       it != estParamsInfo.end(); it++)
-    delete it->prior;
-
-  return logPD;
-}
-
-void
-mexFunction(int nlhs, mxArray *plhs[],
-            int nrhs, const mxArray *prhs[])
-{
-  if (nrhs != 7 )
-    DYN_MEX_FUNC_ERR_MSG_TXT("logposterior: exactly 7 input arguments are required.");
-
-  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);
-
-  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];
-
-  const 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]);
-
-  // Compute and return the value
-  try
-    {
-      *lik = logposterior(estParams, data, options_, M_, estim_params_, bayestopt_, oo_,
-				steadyState, trend_coeff, deepParams, H, Q);
-      *info_mx = 0;
-      *exit_flag = 0;
-    }
-  catch (LogposteriorMexErrMsgTxtException e)
-    {
-      DYN_MEX_FUNC_ERR_MSG_TXT(e.getErrMsg());
-    }
-  catch (SteadyStateSolver::SteadyStateException e)
-    {
-      DYN_MEX_FUNC_ERR_MSG_TXT(e.message.c_str());
-    }
-}

mex/sources/estimation/tests/DsgeLikelihood.m

-function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data,data_index,number_of_observations,no_more_missing_observations)
-% function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data,data_index,number_of_observations,no_more_missing_observations)
-% Evaluates the posterior kernel of a dsge model. 
-% 
-% INPUTS 
-%   xparam1                        [double]   vector of model parameters.
-%   gend                           [integer]  scalar specifying the number of observations.
-%   data                           [double]   matrix of data
-%   data_index                     [cell]     cell of column vectors
-%   number_of_observations         [integer]
-%   no_more_missing_observations   [integer] 
-% OUTPUTS 
-%   fval        :     value of the posterior kernel at xparam1.
-%   cost_flag   :     zero if the function returns a penalty, one otherwise.
-%   ys          :     steady state of original endogenous variables
-%   trend_coeff :
-%   info        :     vector of informations about the penalty:
-%                     41: one (many) parameter(s) do(es) not satisfied the lower bound
-%                     42: one (many) parameter(s) do(es) not satisfied the upper bound
-%               
-% SPECIAL REQUIREMENTS
-%
-
-% Copyright (C) 2004-2012 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/>.
-
-global bayestopt_ estim_params_ options_ trend_coeff_ M_ oo_
-fval            = [];
-ys              = [];
-trend_coeff     = [];
-cost_flag       = 1;
-nobs            = size(options_.varobs,1);
-%------------------------------------------------------------------------------
-% 1. Get the structural parameters & define penalties
-%------------------------------------------------------------------------------
-if options_.mode_compute ~= 1 & any(xparam1 < bayestopt_.lb)
-    k = find(xparam1 < bayestopt_.lb);
-    fval = bayestopt_.penalty+sum((bayestopt_.lb(k)-xparam1(k)).^2);
-    cost_flag = 0;
-    info = 41;
-    return;
-end
-if options_.mode_compute ~= 1 & any(xparam1 > bayestopt_.ub)
-    k = find(xparam1 > bayestopt_.ub);
-    fval = bayestopt_.penalty+sum((xparam1(k)-bayestopt_.ub(k)).^2);
-    cost_flag = 0;
-    info = 42;
-    return;
-end
-Q = M_.Sigma_e;
-H = M_.H;
-for i=1:estim_params_.nvx
-    k =estim_params_.var_exo(i,1);
-    Q(k,k) = xparam1(i)*xparam1(i);
-end
-offset = estim_params_.nvx;
-if estim_params_.nvn
-    for i=1:estim_params_.nvn
-        k = estim_params_.var_endo(i,1);
-        H(k,k) = xparam1(i+offset)*xparam1(i+offset);
-    end
-    offset = offset+estim_params_.nvn;
-end
-if estim_params_.ncx
-    for i=1:estim_params_.ncx
-        k1 =estim_params_.corrx(i,1);
-        k2 =estim_params_.corrx(i,2);
-        Q(k1,k2) = xparam1(i+offset)*sqrt(Q(k1,k1)*Q(k2,k2));
-        Q(k2,k1) = Q(k1,k2);
-    end
-    [CholQ,testQ] = chol(Q);
-    if testQ    %% 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.
-        a = diag(eig(Q));
-        k = find(a < 0);
-        if k > 0
-            fval = bayestopt_.penalty+sum(-a(k));
-            cost_flag = 0;
-            info = 43;
-            return
-        end
-    end
-    offset = offset+estim_params_.ncx;
-end
-if estim_params_.ncn 
-    for i=1:estim_params_.ncn
-        k1 = options_.lgyidx2varobs(estim_params_.corrn(i,1));
-        k2 = options_.lgyidx2varobs(estim_params_.corrn(i,2));
-        H(k1,k2) = xparam1(i+offset)*sqrt(H(k1,k1)*H(k2,k2));
-        H(k2,k1) = H(k1,k2);
-    end
-    [CholH,testH] = chol(H);
-    if testH
-        a = diag(eig(H));
-        k = find(a < 0);
-        if k > 0
-            fval = bayestopt_.penalty+sum(-a(k));
-            cost_flag = 0;
-            info = 44;
-            return
-        end
-    end
-    offset = offset+estim_params_.ncn;
-end
-if estim_params_.np > 0
-    M_.params(estim_params_.param_vals(:,1)) = xparam1(offset+1:end);
-end
-M_.Sigma_e = Q;
-M_.H = H;
-%------------------------------------------------------------------------------
-% 2. call model setup & reduction program
-%------------------------------------------------------------------------------
-[T,R,SteadyState,info] = dynare_resolve(bayestopt_.restrict_var_list,...
-                                        bayestopt_.restrict_columns,...
-                                        bayestopt_.restrict_aux);
-if info(1) == 1 || info(1) == 2 || info(1) == 5
-    fval = bayestopt_.penalty+1;
-    cost_flag = 0;
-    return
-elseif info(1) == 3 || info(1) == 4 || info(1)==6 ||info(1) == 19 || info(1) == 20 || info(1) == 21
-    fval = bayestopt_.penalty+info(2);
-    cost_flag = 0;
-    return
-end
-bayestopt_.mf = bayestopt_.mf1;
-if options_.noconstant
-    constant = zeros(nobs,1);  
-else    
-    if options_.loglinear
-        constant = log(SteadyState(bayestopt_.mfys));
-    else
-        constant = SteadyState(bayestopt_.mfys);
-    end
-end
-if bayestopt_.with_trend
-    trend_coeff = zeros(nobs,1);
-    t = options_.trend_coeffs;
-    for i=1:length(t)
-        if ~isempty(t{i})
-            trend_coeff(i) = evalin('base',t{i});
-        end
-    end
-    trend = repmat(constant,1,gend)+trend_coeff*[1:gend];
-else
-    trend = repmat(constant,1,gend);
-end
-start = options_.presample+1;
-np    = size(T,1);
-mf    = bayestopt_.mf;
-no_missing_data_flag = (number_of_observations==gend*nobs);
-%------------------------------------------------------------------------------
-% 3. Initial condition of the Kalman filter
-%------------------------------------------------------------------------------
-T
-R
-Q
-R*Q*R'
-pause
-options_.lik_init = 1;
-kalman_algo = options_.kalman_algo;
-if options_.lik_init == 1               % Kalman filter
-    if kalman_algo ~= 2
-        kalman_algo = 1;
-    end
-    Pstar = lyapunov_symm(T,R*Q*R',options_.qz_criterium,options_.lyapunov_complex_threshold);
-    Pinf        = [];
-elseif options_.lik_init == 2   % Old Diffuse Kalman filter
-    if kalman_algo ~= 2
-        kalman_algo = 1;
-    end
-    Pstar = options_.Harvey_scale_factor*eye(np);
-    Pinf = [];
-elseif options_.lik_init == 3   % Diffuse Kalman filter
-    if kalman_algo ~= 4
-        kalman_algo = 3;
-    end
-    [QT,ST] = schur(T);
-    e1 = abs(ordeig(ST)) > 2-options_.qz_criterium;
-    [QT,ST] = ordschur(QT,ST,e1);
-    k = find(abs(ordeig(ST)) > 2-options_.qz_criterium);
-    nk = length(k);
-    nk1 = nk+1;
-    Pinf = zeros(np,np);
-    Pinf(1:nk,1:nk) = eye(nk);
-    Pstar = zeros(np,np);
-    B = QT'*R*Q*R'*QT;
-    for i=np:-1:nk+2
-        if ST(i,i-1) == 0
-            if i == np
-                c = zeros(np-nk,1);
-            else
-                c = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i,i+1:end)')+...
-                    ST(i,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i);
-            end
-            q = eye(i-nk)-ST(nk1:i,nk1:i)*ST(i,i);
-            Pstar(nk1:i,i) = q\(B(nk1:i,i)+c);
-            Pstar(i,nk1:i-1) = Pstar(nk1:i-1,i)';
-        else
-            if i == np
-                c = zeros(np-nk,1);
-                c1 = zeros(np-nk,1);
-            else
-                c = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i,i+1:end)')+...
-                    ST(i,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i)+...
-                    ST(i,i-1)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i-1);
-                c1 = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i-1,i+1:end)')+...
-                     ST(i-1,i-1)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i-1)+...
-                     ST(i-1,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i);
-            end
-            q = [eye(i-nk)-ST(nk1:i,nk1:i)*ST(i,i) -ST(nk1:i,nk1:i)*ST(i,i-1);...
-                 -ST(nk1:i,nk1:i)*ST(i-1,i) eye(i-nk)-ST(nk1:i,nk1:i)*ST(i-1,i-1)];
-            z =  q\[B(nk1:i,i)+c;B(nk1:i,i-1)+c1];
-            Pstar(nk1:i,i) = z(1:(i-nk));
-            Pstar(nk1:i,i-1) = z(i-nk+1:end);
-            Pstar(i,nk1:i-1) = Pstar(nk1:i-1,i)';
-            Pstar(i-1,nk1:i-2) = Pstar(nk1:i-2,i-1)';
-            i = i - 1;
-        end
-    end
-    if i == nk+2
-        c = ST(nk+1,:)*(Pstar(:,nk+2:end)*ST(nk1,nk+2:end)')+ST(nk1,nk1)*ST(nk1,nk+2:end)*Pstar(nk+2:end,nk1);
-        Pstar(nk1,nk1)=(B(nk1,nk1)+c)/(1-ST(nk1,nk1)*ST(nk1,nk1));
-    end
-    Z = QT(mf,:);
-    R1 = QT'*R;
-    [QQ,RR,EE] = qr(Z*ST(:,1:nk),0);
-    k = find(abs(diag([RR; zeros(nk-size(Z,1),size(RR,2))])) < 1e-8);
-    if length(k) > 0
-        k1 = EE(:,k);
-        dd =ones(nk,1);
-        dd(k1) = zeros(length(k1),1);
-        Pinf(1:nk,1:nk) = diag(dd);
-    end
-end
-if kalman_algo == 2
-end
-kalman_tol = options_.kalman_tol;
-riccati_tol = options_.riccati_tol;
-mf = bayestopt_.mf1;
-Y   = data-trend;
-Pstar
-pause
-%------------------------------------------------------------------------------
-% 4. Likelihood evaluation
-%------------------------------------------------------------------------------
-if (kalman_algo==1)% Multivariate Kalman Filter
-    if no_missing_data_flag
-        LIK = kalman_filter(T,R,Q,H,Pstar,Y,start,mf,kalman_tol,riccati_tol); 
-    else
-        LIK = ...
-            missing_observations_kalman_filter(T,R,Q,H,Pstar,Y,start,mf,kalman_tol,riccati_tol, ...
-                                               data_index,number_of_observations,no_more_missing_observations);
-    end
-    if isinf(LIK)
-        kalman_algo = 2;
-    end
-end
-if (kalman_algo==2)% Univariate Kalman Filter
-    no_correlation_flag = 1;
-    if length(H)==1 & H == 0
-        H = zeros(nobs,1);
-    else
-        if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal...
-            H = diag(H);
-        else
-            no_correlation_flag = 0;
-        end
-    end
-    if no_correlation_flag
-        LIK = univariate_kalman_filter(T,R,Q,H,Pstar,Y,start,mf,kalman_tol,riccati_tol,data_index,number_of_observations,no_more_missing_observations);
-    else
-        LIK = univariate_kalman_filter_corr(T,R,Q,H,Pstar,Y,start,mf,kalman_tol,riccati_tol,data_index,number_of_observations,no_more_missing_observations);
-    end
-end
-if (kalman_algo==3)% Multivariate Diffuse Kalman Filter
-    if no_missing_data_flag
-        LIK = diffuse_kalman_filter(ST,R1,Q,H,Pinf,Pstar,Y,start,Z,kalman_tol, ...
-                                    riccati_tol);
-    else
-        LIK = missing_observations_diffuse_kalman_filter(ST,R1,Q,H,Pinf, ...
-                                                         Pstar,Y,start,Z,kalman_tol,riccati_tol,...
-                                                         data_index,number_of_observations,...
-                                                         no_more_missing_observations);
-    end
-    if isinf(LIK)
-        kalman_algo = 4;
-    end
-end
-if (kalman_algo==4)% Univariate Diffuse Kalman Filter
-    no_correlation_flag = 1;
-    if length(H)==1 & H == 0
-        H = zeros(nobs,1);
-    else
-        if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal...
-            H = diag(H);
-        else
-            no_correlation_flag = 0;
-        end
-    end
-    if no_correlation_flag
-        LIK = univariate_diffuse_kalman_filter(ST,R1,Q,H,Pinf,Pstar,Y, ...
-                                               start,Z,kalman_tol,riccati_tol,data_index,...
-                                               number_of_observations,no_more_missing_observations);
-    else
-        LIK = univariate_diffuse_kalman_filter_corr(ST,R1,Q,H,Pinf,Pstar, ...
-                                                    Y,start,Z,kalman_tol,riccati_tol,...
-                                                    data_index,number_of_observations,...
-                                                    no_more_missing_observations);
-    end
-end
-if isnan(LIK)
-    cost_flag = 0;
-    return
-end
-if imag(LIK)~=0
-    likelihood = bayestopt_.penalty;
-else
-    likelihood = LIK;
-end
-% ------------------------------------------------------------------------------
-% Adds prior if necessary
-% ------------------------------------------------------------------------------
-lnprior = priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7,bayestopt_.p3,bayestopt_.p4);
-fval    = (likelihood-lnprior);
-likelihood
-lnprior
-fval
-pause
-LIKDLL=logposterior(xparam1,Y,mexext)
-pause
-options_.kalman_algo = kalman_algo;

mex/sources/estimation/tests/Makefile.am

-check_PROGRAMS = test-dr testModelSolution testInitKalman testKalman testPDF
-
-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../../
-
-testModelSolution_SOURCES = ../libmat/Matrix.cc ../libmat/Vector.cc ../libmat/QRDecomposition.cc ../libmat/GeneralizedSchurDecomposition.cc ../libmat/LUSolver.cc ../utils/dynamic_dll.cc ../DecisionRules.cc ../ModelSolution.cc testModelSolution.cc
-testModelSolution_LDADD = $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS) $(LIBADD_DLOPEN)
-testModelSolution_CPPFLAGS = -I.. -I../libmat -I../../ -I../utils
-
-testInitKalman_SOURCES = ../libmat/Matrix.cc ../libmat/Vector.cc ../libmat/QRDecomposition.cc ../libmat/GeneralizedSchurDecomposition.cc ../libmat/LUSolver.cc ../utils/dynamic_dll.cc ../DecisionRules.cc ../ModelSolution.cc ../InitializeKalmanFilter.cc ../DetrendData.cc testInitKalman.cc
-testInitKalman_LDADD = $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS) $(LIBADD_DLOPEN)
-testInitKalman_CPPFLAGS = -I.. -I../libmat -I../../ -I../utils
-
-testKalman_SOURCES = ../libmat/Matrix.cc ../libmat/Vector.cc ../libmat/QRDecomposition.cc ../libmat/GeneralizedSchurDecomposition.cc ../libmat/LUSolver.cc ../utils/dynamic_dll.cc ../DecisionRules.cc ../ModelSolution.cc ../InitializeKalmanFilter.cc ../DetrendData.cc ../KalmanFilter.cc testKalman.cc
-testKalman_LDADD = $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS) $(LIBADD_DLOPEN)
-testKalman_CPPFLAGS = -I.. -I../libmat -I../../ -I../utils
-
-testPDF_SOURCES = ../Prior.cc ../Prior.hh testPDF.cc
-testPDF_CPPFLAGS = -I..
-
-check-local:
-	./test-dr
-	./testPDF

mex/sources/estimation/tests/logposterior_dll_test/random_walk_metropolis_hastings_core.m

-function myoutput = random_walk_metropolis_hastings_core(myinputs,fblck,nblck,whoiam, ThisMatlab)
-% PARALLEL CONTEXT
-% This function contain the most computationally intensive portion of code in
-% random_walk_metropolis_hastings (the 'for xxx = fblck:nblck' loop). The branches in 'for'
-% cycle and are completely independent than suitable to be executed in parallel way.
-% 
-% INPUTS 
-%   o myimput            [struc]     The mandatory variables for local/remote
-%                                    parallel computing obtained from random_walk_metropolis_hastings.m
-%                                    function.
-%   o fblck and nblck    [integer]   The Metropolis-Hastings chains.
-%   o whoiam             [integer]   In concurrent programming a modality to refer to the differents thread running in parallel is needed.
-%                                    The integer whoaim is the integer that
-%                                    allows us to distinguish between them. Then it is the index number of this CPU among all CPUs in the
-%                                    cluster.
-%   o ThisMatlab         [integer]   Allows us to distinguish between the
-%                                    'main' matlab, the slave matlab worker, local matlab, remote matlab,
-%                                     ... Then it is the index number of this slave machine in the cluster.
-% OUTPUTS
-%   o myoutput  [struc]
-%               If executed without parallel is the original output of 'for b =
-%               fblck:nblck' otherwise a portion of it computed on a specific core or
-%               remote machine. In this case:
-%                               record;
-%                               irun;
-%                               NewFile;
-%                               OutputFileName
-%
-% ALGORITHM 
-%   Portion of Metropolis-Hastings.       
-%
-% SPECIAL REQUIREMENTS.
-%   None.
-
-% PARALLEL CONTEXT
-% The most computationally intensive part of this function may be executed
-% in parallel. The code sutable to be executed in parallel on multi core or cluster machine,
-% is removed from this function and placed in random_walk_metropolis_hastings_core.m funtion.
-% Then the DYNARE parallel package contain a set of pairs matlab functios that can be executed in
-% parallel and called name_function.m and name_function_core.m.
-% In addition in the parallel package we have second set of functions used
-% to manage the parallel computation.
-%
-% This function was the first function to be parallelized, later other
-% functions have been parallelized using the same methodology.
-% Then the comments write here can be used for all the other pairs of
-% parallel functions and also for management funtions.
-
-
-% Copyright (C) 2006-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/>.
-
-if nargin<4,
-    whoiam=0;
-end
-
-
-global bayestopt_ estim_params_ options_  M_ oo_
-
-% reshape 'myinputs' for local computation.
-% In order to avoid confusion in the name space, the instruction struct2local(myinputs) is replaced by:
-
-TargetFun=myinputs.TargetFun;
-ProposalFun=myinputs.ProposalFun;
-xparam1=myinputs.xparam1;
-vv=myinputs.vv;
-mh_bounds=myinputs.mh_bounds;
-ix2=myinputs.ix2;
-ilogpo2=myinputs.ilogpo2;
-ModelName=myinputs.ModelName;
-fline=myinputs.fline;
-npar=myinputs.npar;
-nruns=myinputs.nruns;
-NewFile=myinputs.NewFile;
-MAX_nruns=myinputs.MAX_nruns;
-d=myinputs.d;
-InitSizeArray=myinputs.InitSizeArray;                    
-record=myinputs.record;
-varargin=myinputs.varargin;
-
-% Necessary only for remote computing!
-if whoiam
- Parallel=myinputs.Parallel;
- % initialize persistent variables in priordens()
- priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7, ...
-     bayestopt_.p3,bayestopt_.p4,1);
-end
-
-% (re)Set the penalty
-bayestopt_.penalty = Inf;
-
-MhDirectoryName = CheckPath('metropolis');
-
-options_.lik_algo = 1;
-OpenOldFile = ones(nblck,1);
-if strcmpi(ProposalFun,'rand_multivariate_normal')
-    n = npar;
-elseif strcmpi(ProposalFun,'rand_multivariate_student')
-    n = options_.student_degrees_of_freedom;
-end
-% load([MhDirectoryName '/' ModelName '_mh_history.mat'],'record');
-%%%%
-%%%% NOW i run the (nblck-fblck+1) metropolis-hastings chains
-%%%%
-
-
-if any(isnan(bayestopt_.jscale))
-    if exist([ModelName '_optimal_mh_scale_parameter.mat'])% This file is created by mode_compute=6.
-        load([ModelName '_optimal_mh_scale_parameter'])
-        proposal_covariance = d*Scale;
-    else
-        error('mh:: Something is wrong. I can''t figure out the value of the scale parameter.')
-    end
-else
-    proposal_covariance = d*diag(bayestopt_.jscale);
-end
-
-
-jloop=0;
-
-for b = fblck:nblck,
-    jloop=jloop+1;
-    randn('state',record.Seeds(b).Normal);
-    rand('state',record.Seeds(b).Unifor);
-    if (options_.load_mh_file~=0)  & (fline(b)>1) & OpenOldFile(b)
-        load(['./' MhDirectoryName '/' ModelName '_mh' int2str(NewFile(b)) ...
-              '_blck' int2str(b) '.mat'])
-        x2 = [x2;zeros(InitSizeArray(b)-fline(b)+1,npar)];
-        logpo2 = [logpo2;zeros(InitSizeArray(b)-fline(b)+1,1)];
-        OpenOldFile(b) = 0;
-    else
-        x2 = zeros(InitSizeArray(b),npar);
-        logpo2 = zeros(InitSizeArray(b),1);
-    end
-    if exist('OCTAVE_VERSION') || options_.console_mode
-        diary off
-        disp(' ')
-    elseif whoiam
-        %       keyboard;
-        waitbarString = ['Please wait... Metropolis-Hastings (' int2str(b) '/' int2str(options_.mh_nblck) ')...'];
-        %       waitbarTitle=['Metropolis-Hastings ',options_.parallel(ThisMatlab).ComputerName];
-        if options_.parallel(ThisMatlab).Local,
-            waitbarTitle=['Local '];
-        else
-            waitbarTitle=[options_.parallel(ThisMatlab).ComputerName];
-        end        
-        fMessageStatus(0,whoiam,waitbarString, waitbarTitle, options_.parallel(ThisMatlab));
-    else,
-        hh = waitbar(0,['Please wait... Metropolis-Hastings (' int2str(b) '/' int2str(options_.mh_nblck) ')...']);
-        set(hh,'Name','Metropolis-Hastings');
-        
-    end
-    isux = 0;
-    jsux = 0;
-    irun = fline(b);
-    j = 1;
-    while j <= nruns(b)
-        par = feval(ProposalFun, ix2(b,:), proposal_covariance, n);
-        if all( par(:) > mh_bounds(:,1) ) & all( par(:) < mh_bounds(:,2) )
-            try
-                logpost = - feval(TargetFun, par(:),varargin{:});               
-            catch
-                logpost = -inf;
-            end
-            % testing logposterior DLL
-            [junk,logpost1] = logposterior(par(:),varargin{2},mexext);
-            if abs(logpost+logpost1) > 1e-10;
-                disp ([logpost -logpost1])
-            end
-        else
-            logpost = -inf;
-        end
-        if (logpost > -inf) && (log(rand) < logpost-ilogpo2(b))
-            x2(irun,:) = par;
-            ix2(b,:) = par;
-            logpo2(irun) = logpost; 
-            ilogpo2(b) = logpost;
-            isux = isux + 1;
-            jsux = jsux + 1;
-        else    
-            x2(irun,:) = ix2(b,:);
-            logpo2(irun) = ilogpo2(b);
-        end
-        prtfrc = j/nruns(b);
-        if exist('OCTAVE_VERSION') || options_.console_mode
-            if mod(j, 10) == 0
-                if exist('OCTAVE_VERSION')
-                    printf('MH: Computing Metropolis-Hastings (chain %d/%d): %3.f%% done, acception rate: %3.f%%\r', b, nblck, 100 * prtfrc, 100 * isux / j);
-                else
-                    fprintf('   MH: Computing Metropolis-Hastings (chain %d/%d): %3.f \b%% done, acceptance rate: %3.f \b%%\r', b, nblck, 100 * prtfrc, 100 * isux / j);
-                end
-            end
-            if mod(j,50)==0 & whoiam  
-                %             keyboard;
-                waitbarString = [ '(' int2str(b) '/' int2str(options_.mh_nblck) '), ' sprintf('accept. %3.f%%%%', 100 * isux/j)];
-                fMessageStatus(prtfrc,whoiam,waitbarString, '', options_.parallel(ThisMatlab));
-            end
-        else
-            if mod(j, 3)==0 & ~whoiam
-                waitbar(prtfrc,hh,[ '(' int2str(b) '/' int2str(options_.mh_nblck) ') ' sprintf('%f done, acceptation rate %f',prtfrc,isux/j)]);
-            elseif mod(j,50)==0 & whoiam,  
-                %             keyboard;
-                waitbarString = [ '(' int2str(b) '/' int2str(options_.mh_nblck) ') ' sprintf('%f done, acceptation rate %f',prtfrc,isux/j)];
-                fMessageStatus(prtfrc,whoiam,waitbarString, waitbarTitle, options_.parallel(ThisMatlab));
-            end
-        end
-        
-        if (irun == InitSizeArray(b)) | (j == nruns(b)) % Now I save the simulations
-            save([MhDirectoryName '/' ModelName '_mh' int2str(NewFile(b)) '_blck' int2str(b) '.mat'],'x2','logpo2');
-            fidlog = fopen([MhDirectoryName '/metropolis.log'],'a');
-            fprintf(fidlog,['\n']);
-            fprintf(fidlog,['%% Mh' int2str(NewFile(b)) 'Blck' int2str(b) ' (' datestr(now,0) ')\n']);
-            fprintf(fidlog,' \n');
-            fprintf(fidlog,['  Number of simulations.: ' int2str(length(logpo2)) '\n']);
-            fprintf(fidlog,['  Acceptation rate......: ' num2str(jsux/length(logpo2)) '\n']);
-            fprintf(fidlog,['  Posterior mean........:\n']);
-            for i=1:length(x2(1,:))
-                fprintf(fidlog,['    params:' int2str(i) ': ' num2str(mean(x2(:,i))) '\n']);
-            end
-            fprintf(fidlog,['    log2po:' num2str(mean(logpo2)) '\n']);
-            fprintf(fidlog,['  Minimum value.........:\n']);;
-            for i=1:length(x2(1,:))
-                fprintf(fidlog,['    params:' int2str(i) ': ' num2str(min(x2(:,i))) '\n']);
-            end
-            fprintf(fidlog,['    log2po:' num2str(min(logpo2)) '\n']);
-            fprintf(fidlog,['  Maximum value.........:\n']);
-            for i=1:length(x2(1,:))
-                fprintf(fidlog,['    params:' int2str(i) ': ' num2str(max(x2(:,i))) '\n']);
-            end
-            fprintf(fidlog,['    log2po:' num2str(max(logpo2)) '\n']);
-            fprintf(fidlog,' \n');
-            fclose(fidlog);
-            jsux = 0;
-            if j == nruns(b) % I record the last draw...
-                record.LastParameters(b,:) = x2(end,:);
-                record.LastLogLiK(b) = logpo2(end);
-            end
-            % size of next file in chain b
-            InitSizeArray(b) = min(nruns(b)-j,MAX_nruns);
-            % initialization of next file if necessary
-            if InitSizeArray(b)
-                x2 = zeros(InitSizeArray(b),npar);
-                logpo2 = zeros(InitSizeArray(b),1);
-                NewFile(b) = NewFile(b) + 1;
-                irun = 0;
-            end
-        end
-        j=j+1;
-        irun = irun + 1;
-    end% End of the simulations for one mh-block.
-    record.AcceptationRates(b) = isux/j;
-    if exist('OCTAVE_VERSION') || options_.console_mode
-        if exist('OCTAVE_VERSION')
-            printf('\n');
-        else
-            fprintf('\n');
-        end
-        diary on;
-    elseif ~whoiam
-        close(hh);
-    end
-    record.Seeds(b).Normal = randn('state');
-    record.Seeds(b).Unifor = rand('state');
-    OutputFileName(jloop,:) = {[MhDirectoryName,filesep], [ModelName '_mh*_blck' int2str(b) '.mat']};
-end% End of the loop over the mh-blocks.
-
-
-myoutput.record = record;
-myoutput.irun = irun;
-myoutput.NewFile = NewFile;
-myoutput.OutputFileName = OutputFileName;
\ No newline at end of file

mex/sources/estimation/tests/logposterior_dll_test/rawdata_euromodel_1.m

-C =[
-   -7.4073
-   -6.1860
-   -6.5983
-   -5.6088
-   -5.0547
-   -4.4774
-   -3.8081
-   -3.8425
-   -2.4178
-   -1.9835
-   -1.0395
-   -0.1583
-   -0.0397
-    0.3505
-   -0.1879
-   -0.0067
-    0.0478
-   -1.2247
-   -1.4349
-   -0.7973
-   -0.0461
-    0.5844
-    1.1372
-    1.3801
-    1.8023
-    2.2972
-    2.0469
-    2.5435
-    2.8169
-    3.2007
-    2.6705
-    3.0518
-    3.2445
-    3.8443
-    3.8525
-    4.9494
-    4.2770
-    4.9532
-    5.1441
-    3.7124
-    3.9880
-    3.6926
-    2.6005
-    1.8679
-    1.9085
-    1.5563
-    1.2308
-    0.3264
-   -0.2208
-   -0.2483
-   -0.4082
-   -1.0315
-   -1.6030
-   -1.5499
-   -1.3777
-   -2.1675
-   -2.5138
-   -2.8820
-   -2.6958
-   -2.4719
-   -1.9854
-   -1.7954
-   -2.2362
-   -1.0595
-   -0.8808
-   -0.8548
-   -1.2839
-   -0.1363
-    0.2104
-    0.8810
-    0.3555
-    0.4766
-    1.3269
-    1.4506
-    1.4308
-    1.6263
-    1.9842
-    2.3948
-    2.8710
-    3.0177
-    2.9305
-    3.1739
-    3.7380
-    3.8285
-    3.3342
-    3.7447
-    3.7830
-    3.1039
-    2.8413
-    3.0338
-    0.3669
-    0.0847
-    0.0104
-    0.2115
-   -0.6649
-   -0.9625
-   -0.7330
-   -0.8664
-   -1.4441
-   -1.0179
-   -1.2729
-   -1.9539
-   -1.4427
-   -2.0371
-   -1.9764
-   -2.5654
-   -2.8570
-   -2.5842
-   -3.0427
-   -2.8312
-   -2.3320
-   -2.2768
-   -2.1816
-   -2.1043
-   -1.8969
-   -2.2388
-   -2.1679
-   -2.1172
-];
-
-E =[
-    0.6263
-    0.7368
-    0.7477
-    1.0150
-    0.6934
-    0.4135
-    0.3845
-    0.2380
-    0.2853
-    0.5999
-    0.8622
-    1.2116
-    1.4921
-    1.5816
-    1.7259
-    1.6276
-    1.2422
-    0.8084
-    0.4710
-   -0.3704
-   -0.6427
-   -0.5323
-   -0.5562
-   -0.3651
-   -0.4356
-   -0.7164
-   -0.5816
-   -0.4635
-   -0.8456
-   -0.9708
-   -0.7138
-   -0.7499
-   -0.6941
-   -0.6656
-   -0.2912
-   -0.1650
-    0.0774
-    0.2307
-    0.4484
-    0.4942
-    0.4653
-    0.2196
-    0.1736
-   -0.1595
-   -0.3918
-   -0.4611
-   -0.8493
-   -0.7384
-   -1.0604
-   -1.2166
-   -1.7187
-   -1.6932
-   -1.7830
-   -1.7035
-   -2.2079
-   -2.3769
-   -2.2511
-   -2.1093
-   -2.4638
-   -2.4027
-   -2.1313
-   -1.9199
-   -1.7941
-   -1.4661
-   -1.2269
-   -1.0392
-   -1.0725
-   -0.7156
-   -0.4778
-   -0.4233
-   -0.0409
-    0.1620
-    0.4280
-    0.5873
-    1.0323
-    1.3420
-    1.6902
-    2.0680
-    2.8219
-    3.2511
-    3.2930
-    3.5633
-    3.8992
-    3.6874
-    3.2849
-    3.1614
-    2.6221
-    2.5067
-    1.9223
-    1.1777
-    0.4483
-   -0.0661
-   -0.4424
-   -0.9000
-   -1.1478
-   -1.2047
-   -1.1412
-   -1.2383
-   -1.1048
-   -0.9716
-   -0.9287
-   -1.0057
-   -1.0827
-   -1.0200
-   -1.0072
-   -1.1740
-   -1.2809
-   -1.1086
-   -0.9866
-   -0.8947
-   -0.5875
-   -0.2329
-    0.1493
-    0.4906
-    0.8400
-    1.0720
-    1.2648
-    1.5431
-];
-
-I =[
-    2.6617
-    2.4325
-    1.9592
-    3.2530
-    2.9949
-    3.7918
-    4.7444
-    4.8289
-    5.5983
-    7.8923
-    9.4297
-    9.5010
-   10.0150
-   10.0413
-    9.6046
-    6.4766
-    5.9647
-    3.0114
-    0.5683
-   -2.1226
-   -2.1855
-   -0.8329
-   -1.5207
-   -1.3419
-   -1.7897
-   -0.1476
-    0.4675
-   -1.6516
-   -1.5419
-   -1.3050
-   -1.2451
-   -0.7815
-   -0.7796
-   -0.3612
-   -2.4072
-    1.1162
-    1.1383
-    3.4132
-    5.0356
-    2.8016
-    2.1734
-    0.9366
-   -0.7050
-   -1.5021
-   -2.9868
-   -6.0237
-   -6.2589
-   -6.9138
-   -8.2340
-   -9.2589
-   -9.2465
-   -9.6988
-   -9.7782
-  -10.5645
-  -10.7544
-  -13.1583
-  -12.2718
-  -12.0131
-  -13.5983
-  -12.3579
-  -10.9146
-  -11.1572
-  -12.4935
-   -9.4393
-   -8.5535
-   -7.3723
-  -10.0169
-   -6.6088
-   -5.2045
-   -4.1024
-   -2.8472
-   -1.3139
-    0.0477
-    1.5629
-    3.6947
-    4.0327
-    4.1320
-    7.1400
-    9.1036
-    8.5609
-    7.6576
-    8.8022
-    8.9611
-   10.0871
-    9.4797
-    9.3964
-   10.0363
-    8.6340
-    6.6522
-    4.4471
-    0.2854
-   -2.1879
-   -2.9879
-   -4.1021
-   -2.7713
-   -2.2281
-   -1.2908
-   -0.3250
-    0.6534
-    0.3942
-    0.3534
-   -0.1532
-   -1.7936
-    0.4909
-    0.3634
-    0.4290
-   -0.9709
-    0.1942
-    0.6103
-    1.4426
-    2.7225
-    1.7525
-    3.2780
-    3.5985
-    4.9011
-    5.3312
-    6.4402
-    6.6529
-];
-
-L =[
-    0.6263
-    0.7368
-    0.7477
-    1.0150
-    0.6934
-    0.4135
-    0.3845
-    0.2380
-    0.2853
-    0.5999
-    0.8622
-    1.2116
-    1.4921
-    1.5816
-    1.7259
-    1.6276
-    1.2422
-    0.8084
-    0.4710
-   -0.3704
-   -0.6427
-   -0.5323
-   -0.5562
-   -0.3651
-   -0.4356
-   -0.7164
-   -0.5816
-   -0.4635
-   -0.8456
-   -0.9708
-   -0.7138
-   -0.7499
-   -0.6941
-   -0.6656
-   -0.2912
-   -0.1650
-    0.0774
-    0.2307
-    0.4484
-    0.4942
-    0.4653
-    0.2196
-    0.1736
-   -0.1595
-   -0.3918
-   -0.4611
-   -0.8493
-   -0.7384
-   -1.0604
-   -1.2166
-   -1.7187
-   -1.6932
-   -1.7830
-   -1.7035
-   -2.2079
-   -2.3769
-   -2.2511
-   -2.1093
-   -2.4638
-   -2.4027
-   -2.1313
-   -1.9199
-   -1.7941
-   -1.4661
-   -1.2269
-   -1.0392
-   -1.0725
-   -0.7156
-   -0.4778
-   -0.4233
-   -0.0409
-    0.1620
-    0.4280
-    0.5873
-    1.0323
-    1.3420
-    1.6902
-    2.0680
-    2.8219
-    3.2511
-    3.2930
-    3.5633
-    3.8992
-    3.6874
-    3.2849
-    3.1614
-    2.6221
-    2.5067
-    1.9223
-    1.1777
-    0.4483
-   -0.0661
-   -0.4424
-   -0.9000
-   -1.1478
-   -1.2047
-   -1.1412
-   -1.2383
-   -1.1048
-   -0.9716
-   -0.9287
-   -1.0057
-   -1.0827
-   -1.0200
-   -1.0072
-   -1.1740
-   -1.2809
-   -1.1086
-   -0.9866
-   -0.8947
-   -0.5875
-   -0.2329
-    0.1493
-    0.4906
-    0.8400
-    1.0720
-    1.2648
-    1.5431
-];
-
-PIE =[
-   -1.0113
-   -0.8305
-    0.2332
-   -0.8746
-   -0.7978
-   -0.9220
-   -0.2487
-   -0.7749
-   -0.5460
-   -0.5347
-    0.5050
-   -0.0334
-    0.6756
-    0.8791
-    0.7267
-    1.0997
-    1.1750
-    1.1927
-    0.4420
-    0.5357
-    0.0345
-    0.0196
-    0.3371
-    0.9379
-    1.2160
-    0.3393
-    0.5813
-    0.7410
-    0.3374
-    0.2616
-    0.4025
-    0.4799
-    0.5981
-   -0.1523
-    0.4458
-    0.2182
-    0.9793
-    0.7562
-    1.0064
-    0.8203
-    0.6966
-    0.3352
-    0.6581
-    0.6111
-    0.9833
-    1.1991
-    0.9562
-    0.3868
-    0.2939
-    0.2471
-    0.8331
-    0.0715
-    0.3910
-    0.3301
-    0.2547
-   -0.2702
-   -0.2998
-   -0.1953
-   -0.2293
-   -0.3284
-    0.0480
-   -0.0374
-    0.3253
-   -0.3434
-   -0.3892
-   -0.7178
-   -0.4758
-   -0.6794
-   -0.8505
-   -0.3512
-   -0.4436
-   -0.5101
-   -0.4574
-   -0.2696
-   -0.1047
-   -0.5745
-   -0.2989
-   -0.0063
-    0.0088
-   -0.1184
-   -0.1506
-   -0.4073
-    0.2674
-    0.2896
-    0.0669
-    0.1166
-   -0.1699
-   -0.2518
-   -0.0562
-   -0.3269
-   -0.0703
-   -0.1046
-   -0.4888
-   -0.3524
-   -0.2485
-   -0.5870
-   -0.4546
-   -0.3970
-   -0.2353
-   -0.0352
-   -0.2171
-   -0.3754
-   -0.4322
-   -0.4572
-   -0.4903
-   -0.4518
-   -0.6435
-   -0.6304
-   -0.4148
-   -0.2892
-   -0.4318
-   -0.6010
-   -0.4148
-   -0.4315
-   -0.3531
-   -0.8053
-   -0.4680
-   -0.4263
-];
-
-R =[
-   -1.0750
-   -1.1540
-   -1.3682
-   -1.4569
-   -1.3490
-   -1.4011
-   -1.6486
-   -1.6968
-   -1.6976
-   -1.2567
-   -1.1392
-   -0.7783
-   -0.3021
-   -0.0435
-    0.0066
-   -0.0043
-    0.1029
-   -0.0628
-   -0.5358
-   -0.9627
-   -1.1079
-   -1.0918
-   -0.9966
-   -0.6223
-   -0.3616
-   -0.2711
-   -0.0997
-   -0.2810
-   -0.3710
-   -0.3167
-   -0.5301
-   -0.5826
-   -0.3194
-   -0.2713
-   -0.5287
-   -0.2432
-    0.1098
-    0.5349
-    0.7094
-    0.8415
-    0.6226
-    0.7376
-    0.9316
-    1.4370
-    1.5853
-    1.4267
-    1.1783
-    1.2046
-    0.9689
-    0.7918
-    0.6315
-    0.5950
-    0.6853
-    0.7171
-    0.5887
-    0.4873
-    0.4027
-    0.3489
-    0.2934
-    0.3060
-    0.1741
-    0.0348
-    0.0771
-   -0.1005
-   -0.1518
-   -0.1104
-   -0.0681
-   -0.0059
-    0.0256
-    0.0404
-   -0.1721
-   -0.2002
-    0.0015
-    0.1249
-    0.3738
-    0.4320
-    0.5579
-    0.8186
-    0.8727
-    0.7356
-    0.7243
-    0.8635
-    0.9058
-    0.7656
-    0.7936
-    0.8631
-    0.9074
-    0.9547
-    1.2045
-    1.0850
-    0.9178
-    0.5242
-    0.3178
-    0.1472
-    0.0227
-   -0.0799
-   -0.0611
-   -0.0140
-    0.1132
-    0.1774
-    0.0782
-    0.0436
-   -0.1596
-   -0.2691
-   -0.2895
-   -0.3791
-   -0.4020
-   -0.4166
-   -0.4037
-   -0.3636
-   -0.4075
-   -0.4311
-   -0.4470
-   -0.5111
-   -0.6274
-   -0.7261
-   -0.6974
-   -0.5012
-];
-
-W =[
-  -14.8791
-  -13.2300
-  -13.5037
-  -13.0249
-  -11.2546
-  -10.0148
-   -8.8586
-   -8.5739
-   -7.7851
-   -6.7136
-   -5.5878
-   -4.6881
-   -3.8039
-   -3.0366
-   -2.7342
-   -1.3135
-   -0.7387
-   -0.1131
-   -0.2769
-    0.8696
-    1.8855
-    2.3667
-    2.4942
-    3.2049
-    3.9682
-    5.1500
-    4.7047
-    4.7827
-    5.3377
-    5.6614
-    5.2813
-    5.2967
-    5.5175
-    6.1526
-    5.6627
-    6.0694
-    6.5824
-    6.9032
-    6.7849
-    6.6896
-    6.6201
-    6.9933
-    5.8959
-    6.7419
-    6.9999
-    6.4009
-    5.5083
-    5.1054
-    5.2813
-    4.5790
-    3.9589
-    3.8599
-    3.8978
-    2.7957
-    3.2480
-    1.4634
-    1.9219
-    1.8398
-    1.9279
-    1.8316
-    1.6092
-    1.2741
-    0.2031
-   -0.0236
-   -0.1004
-   -0.3034
-   -1.0273
-   -0.2205
-    0.0458
-    0.2386
-   -0.0977
-   -0.3145
-   -0.1416
-   -0.7009
-   -0.9082
-   -0.8802
-   -0.5644
-   -0.5852
-   -0.5346
-    0.0652
-    0.1301
-    0.3444
-   -0.3592
-    0.8096
-    0.9644
-    1.0289
-    1.2781
-    1.2298
-    2.2134
-    2.0808
-    0.4925
-    0.6506
-    0.5531
-    0.2456
-   -0.5351
-   -0.8183
-   -0.8967
-   -0.7268
-   -1.0738
-   -1.2844
-   -1.4338
-   -1.6995
-   -1.7085
-   -2.2889
-   -2.1018
-   -2.4273
-   -2.4609
-   -2.1407
-   -2.3847
-   -3.1689
-   -4.5581
-   -4.1027
-   -4.2436
-   -4.8836
-   -5.9660
-   -4.9971
-   -5.2386
-   -5.6618
-];
-
-Y =[
-   -4.9347
-   -4.6205
-   -5.2198
-   -4.5937
-   -3.8015
-   -3.6643
-   -2.7239
-   -2.7524
-   -2.0634
-   -1.0112
-    0.0530
-    0.7623
-    1.7927
-    2.1486
-    2.4866
-    2.1456
-    2.1671
-   -0.0254
-   -1.6716
-   -1.9673
-   -1.6109
-   -1.0292
-   -0.1222
-    0.7329
-    1.1234
-    2.0603
-    1.7998
-    1.4820
-    1.1732
-    1.6424
-    1.5382
-    2.1399
-    2.0127
-    2.7210
-    2.4966
-    3.5249
-    3.6237
-    4.2011
-    4.5634
-    3.3442
-    2.7761
-    1.9812
-    1.3779
-    1.4616
-    1.3029
-    0.7594
-    0.3695
-    0.0832
-   -0.8118
-   -1.4557
-   -1.4850
-   -1.2346
-   -1.5696
-   -1.3785
-   -0.7682
-   -2.0308
-   -1.7778
-   -1.7801
-   -2.1711
-   -1.7469
-   -1.3413
-   -1.3352
-   -2.4390
-   -1.2125
-   -1.1695
-   -1.0891
-   -2.4753
-   -1.3503
-   -0.9412
-   -0.1470
-    0.0026
-    0.1108
-    0.6890
-    1.3520
-    1.6018
-    2.0667
-    1.7625
-    2.6658
-    3.4048
-    3.2507
-    3.4251
-    3.2174
-    3.1903
-    3.3396
-    3.1358
-    2.8625
-    3.3546
-    2.4609
-    1.9534
-    0.9962
-   -0.7904
-   -1.1672
-   -1.2586
-   -1.3593
-   -1.3443
-   -0.9413
-   -0.6023
-   -0.4516
-   -0.5129
-   -0.8741
-   -1.0784
-   -1.4091
-   -1.3627
-   -1.5731
-   -1.6037
-   -1.8814
-   -2.1482
-   -1.3597
-   -1.1855
-   -1.1122
-   -0.8424
-   -0.9747
-   -1.1385
-   -1.4548
-   -1.4284
-   -1.4633
-   -1.0621
-   -0.7871
-];

mex/sources/estimation/tests/logposterior_dll_test/sweuromodel_dll.mod

-//options_.usePartInfo=1;
-
-var MC E EF R_KF QF CF IF YF LF PIEF WF RF R_K Q C I Y L PIE W R EE_A PIE_BAR EE_B EE_G EE_L EE_I KF K one BIGTHETA;    
-
-varexo E_A E_B E_G E_L E_I ETA_R E_PIE_BAR ETA_Q ETA_P ETA_W  ;  
- 
-parameters 
-xi_e lambda_w alpha czcap beta phi_i tau sig_c hab ccs cinvs phi_y gamma_w xi_w gamma_p xi_p sig_l r_dpi 
-r_pie r_dy r_y rho rho_a rho_pb rho_b rho_g rho_l rho_i LMP  ;
-
-
-
-alpha=.30;
-beta=.99;
-tau=0.025;
-ccs=0.6;
-cinvs=.22;  //% alpha*(tau+ctrend)/R_K   R_K=ctrend/beta-1+tau  
-lambda_w = 0.5;
-phi_i= 6.771;
-sig_c=   1.353; 
-hab=    0.573;    
-xi_w=   0.737;
-sig_l=    2.400;
-xi_p=   0.908;
-xi_e= 0.599;
-gamma_w=    0.763;
-gamma_p=    0.469;
-czcap=    0.169;
-phi_y=    1.408;
-r_pie=     1.684;
-r_dpi=    0.14;
-rho=      0.961;
-r_y=      0.099;
-r_dy=     0.159;
-rho_a=    0.823;
-rho_b=    0.855;
-rho_g=    0.949;
-rho_l=   0.889;
-rho_i=   0.927;
-rho_pb=  0.924;
-LMP = 0.0 ; //NEW.
-
-model(linear, use_dll); 
-          CF = (1/(1+hab))*(CF(1)+hab*CF(-1))-((1-hab)/((1+hab)*sig_c))*(RF-PIEF(1)-EE_B) ;
-	      0 =  alpha*R_KF+(1-alpha)*WF -EE_A ;
-          PIEF = 0*one;
-	      IF = (1/(1+beta))* ((  IF(-1) + beta*(IF(1)))+(1/phi_i)*QF)+0*ETA_Q+EE_I ;
-	      QF = -(RF-PIEF(1))+(1-beta*(1-tau))*((0+czcap)/czcap)*R_KF(1)+beta*(1-tau)*QF(1) +0*EE_I ;
-          KF =  (1-tau)*KF(-1)+tau*IF(-1) ;
-	      YF = (ccs*CF+cinvs*IF)+EE_G  ;
-
-	      YF = 1*phi_y*( alpha*KF+alpha*(1/czcap)*R_KF+(1-alpha)*LF+EE_A ) ;
-	      WF = (sig_c/(1-hab))*(CF-hab*CF(-1)) + sig_l*LF - EE_L ;
-	      LF = R_KF*((1+czcap)/czcap)-WF+KF ;
-          EF = EF(-1)+EF(1)-EF+(LF-EF)*((1-xi_e)*(1-xi_e*beta)/(xi_e));
-         
-	      C = (hab/(1+hab))*C(-1)+(1/(1+hab))*C(1)-((1-hab)/((1+hab)*sig_c))*(R-PIE(1)-EE_B) ;
-	      I = (1/(1+beta))* ((  I(-1) + beta*(I(1)))+(1/phi_i)*Q )+1*ETA_Q+1*EE_I ;
-	      Q = -(R-PIE(1))+(1-beta*(1-tau))*((0+czcap)/czcap)*R_K(1)+beta*(1-tau)*Q(1) +EE_I*0+0*ETA_Q ;
-	      K =  (1-tau)*K(-1)+tau*I(-1) ;
-	      Y = (ccs*C+cinvs*I)+ EE_G   ;
-	      Y = phi_y*( alpha*K+alpha*(1/czcap)*R_K+(1-alpha)*L ) +phi_y*EE_A  ;
-	      PIE = (1/(1+beta*gamma_p))*
-	            ( 
-	            (beta)*(PIE(1)) +(gamma_p)*(PIE(-1)) 
-	            +((1-xi_p)*(1-beta*xi_p)/(xi_p))*(MC)
-	            )  + ETA_P ; 
-	            
-	      MC = alpha*R_K+(1-alpha)*W -EE_A;
-	      W =  (1/(1+beta))*(beta*W(+1)+W(-1))
-                +(beta/(1+beta))*(PIE(+1))
-                -((1+beta*gamma_w)/(1+beta))*(PIE)
-                +(gamma_w/(1+beta))*(PIE(-1))
-                -(1/(1+beta))*(((1-beta*xi_w)*(1-xi_w))/(((1+(((1+lambda_w)*sig_l)/(lambda_w))))*xi_w))*(W-sig_l*L-(sig_c/(1-hab))*(C-hab*C(-1))+EE_L)
-                +ETA_W;
-	      L = R_K*((1+czcap)/czcap)-W+K ;
-
-//	      R = r_dpi*(PIE-PIE(-1))
-//              +(1-rho)*(r_pie*(PIE(-1)-PIE_BAR)+r_y*(Y-YF))
-//              +r_dy*(Y-YF-(Y(-1)-YF(-1)))
-//              +rho*(R(-1)-PIE_BAR)
-//              +PIE_BAR
-//              +ETA_R;
-
-
-	      R = 
-
-r_dpi*(PIE-PIE(-1))
-
-              +(1-rho)*(r_pie*(BIGTHETA)+r_y*(Y-YF))
-              +r_dy*(Y-YF-(Y(-1)-YF(-1)))
-              +rho*(R(-1)-PIE_BAR)
-              +PIE_BAR
-              +ETA_R;
-
-
-          E = E(-1)+E(1)-E+(L-E)*((1-xi_e)*(1-xi_e*beta)/(xi_e));
-          
-          
-          EE_A = (rho_a)*EE_A(-1)  + E_A;
-	      PIE_BAR = rho_pb*PIE_BAR(-1)+ E_PIE_BAR ;
-	      EE_B = rho_b*EE_B(-1) + E_B ;
-	      EE_G = rho_g*EE_G(-1) + E_G ;
-	      EE_L = rho_l*EE_L(-1) + E_L ;
-	      EE_I = rho_i*EE_I(-1) + E_I ;
-	      one = 0*one(-1) ;
-
-		LMP*BIGTHETA(1) = BIGTHETA - (PIE(-1) - PIE_BAR) ; 
-
-end; 
-
- 
-shocks;
-var E_A; stderr 0.598;
-var E_B; stderr 0.336;
-var E_G; stderr 0.325;
-var E_I; stderr 0.085;
-var E_L; stderr 3.520;
-var ETA_P; stderr 0.160;
-var ETA_W; stderr 0.289;
-var ETA_R; stderr 0.081;
-var ETA_Q; stderr 0.604;
-var E_PIE_BAR; stderr 0.017;
-end;
-
-//stoch_simul(irf=20) Y C PIE R W R_K L Q I K ;
-
-// stoch_simul generates what kind of standard errors for the shocks ?
-
-//steady;
-//check;
-//stoch_simul(periods=200,irf=20,simul_seed=3) Y C PIE MC R W R_K E L I ;
-
-//datatomfile('ddd',[]);
-
-// new syntax 
-
-estimated_params;
-// PARAM NAME, INITVAL, LB, UB, PRIOR_SHAPE, PRIOR_P1, PRIOR_P2, PRIOR_P3, PRIOR_P4, JSCALE
-// PRIOR_SHAPE: BETA_PDF, GAMMA_PDF, NORMAL_PDF, INV_GAMMA_PDF
-stderr E_A,0.543,0.01,4,INV_GAMMA_PDF,0.4,2;
-stderr E_PIE_BAR,0.072,0.001,4,INV_GAMMA_PDF,0.02,10;
-stderr E_B,0.2694,0.01,4,INV_GAMMA_PDF,0.2,2;
-stderr E_G,0.3052,0.01,4,INV_GAMMA_PDF,0.3,2;
-stderr E_L,1.4575,0.1,6,INV_GAMMA_PDF,1,2;
-stderr E_I,0.1318,0.01,4,INV_GAMMA_PDF,0.1,2;
-stderr ETA_R,0.1363,0.01,4,INV_GAMMA_PDF,0.1,2;
-stderr ETA_Q,0.4842,0.01,4,INV_GAMMA_PDF,0.4,2;
-stderr ETA_P,0.1731,0.01,4,INV_GAMMA_PDF,0.15,2;
-stderr ETA_W,0.2462,0.1,4,INV_GAMMA_PDF,0.25,2;
-rho_a,.9722,.1,.9999,BETA_PDF,0.85,0.1;
-rho_pb,.85,.1,.999,BETA_PDF,0.85,0.1;
-rho_b,.7647,.1,.99,BETA_PDF,0.85,0.1;
-rho_g,.9502,.1,.9999,BETA_PDF,0.85,0.1;
-rho_l,.9542,.1,.9999,BETA_PDF,0.85,0.1;
-rho_i,.6705,.1,.99,BETA_PDF,0.85,0.1;
-phi_i,5.2083,1,15,NORMAL_PDF,4,1.5;
-sig_c,0.9817,0.25,3,NORMAL_PDF,1,0.375;
-hab,0.5612,0.3,0.95,BETA_PDF,0.7,0.1;
-xi_w,0.7661,0.3,0.9,BETA_PDF,0.75,0.05;
-sig_l,1.7526,0.5,5,NORMAL_PDF,2,0.75;
-xi_p,0.8684,0.3,0.95,BETA_PDF,0.75,0.05;
-xi_e,0.5724,0.1,0.95,BETA_PDF,0.5,0.15;
-gamma_w,0.6202,0.1,0.99,BETA_PDF,0.75,0.15;
-gamma_p,0.6638,0.1,0.99,BETA_PDF,0.75,0.15;
-czcap,0.2516,0.01,2,NORMAL_PDF,0.2,0.075;
-phi_y,1.3011,1.001,2,NORMAL_PDF,1.45,0.125;
-r_pie,1.4616,1.2,2,NORMAL_PDF,1.7,0.1;
-r_dpi,0.1144,0.01,0.5,NORMAL_PDF,0.3,0.1;
-rho,0.8865,0.5,0.99,BETA_PDF,0.8,0.10;
-r_y,0.0571,0.01,0.2,NORMAL_PDF,0.125,0.05;
-r_dy,0.2228,0.05,0.5,NORMAL_PDF,0.0625,0.05;
-end;
-
-varobs Y C I E PIE W R;
-
-//estimation(datafile=rawdata_euromodel_1,presample=40, first_obs=1, nobs=118, lik_init=2, mode_compute=1,mh_replic=0);
-estimation(datafile=rawdata_euromodel_1,presample=40, first_obs=1, nobs=118,mh_jscale=0.2,mh_replic=1000
-//,mode_compute=0,mode_file=sweuromodel_dll_mode
-);
-
-
-//stoch_simul(periods=200,irf=20,simul_seed=3) Y C PIE R W R_K L Q I K ;
-