Added new loglikelihood DLL (does not yet contain prior computation, only the likelihood)

mex/build/loglikelihood.am

+vpath %.cc $(top_srcdir)/../../sources/estimation $(top_srcdir)/../../sources/estimation/libmat $(top_srcdir)/../../sources/estimation/utils
+vpath %.hh $(top_srcdir)/../../sources/estimation $(top_srcdir)/../../sources/estimation/libmat
+
+CPPFLAGS += -I$(top_srcdir)/../../sources/estimation/libmat -I$(top_srcdir)/../../sources/estimation/utils
+
+noinst_PROGRAMS = loglikelihood
+
+loglikelihood_LDADD = $(LIBADD_DLOPEN)
+
+MAT_SRCS = \
+	Matrix.hh \
+	Matrix.cc \
+	Vector.hh \
+	Vector.cc \
+	BlasBindings.hh \
+	DiscLyapFast.hh \
+	GeneralizedSchurDecomposition.cc \
+	GeneralizedSchurDecomposition.hh \
+	LapackBindings.hh \
+	LUSolver.cc \
+	LUSolver.hh \
+	QRDecomposition.cc \
+	QRDecomposition.hh \
+	VDVEigDecomposition.cc \
+	VDVEigDecomposition.hh
+
+nodist_loglikelihood_SOURCES = \
+	$(MAT_SRCS) \
+	DecisionRules.cc \
+	DecisionRules.hh \
+	DetrendData.cc \
+	DetrendData.hh \
+	EstimatedParameter.cc \
+	EstimatedParameter.hh \
+	EstimatedParametersDescription.cc \
+	EstimatedParametersDescription.hh \
+	EstimationSubsample.cc \
+	EstimationSubsample.hh \
+	InitializeKalmanFilter.cc \
+	InitializeKalmanFilter.hh \
+	KalmanFilter.cc \
+	KalmanFilter.hh \
+	LogLikelihoodSubSample.cc \
+	LogLikelihoodSubSample.hh \
+	LogLikelihoodMain.hh \
+	LogLikelihoodMain.cc \
+	ModelSolution.cc \
+	ModelSolution.hh \
+	Prior.cc \
+	Prior.hh \
+	dynamic_dll.cc \
+	dynamic_dll.hh \
+	loglikelihood.cc

mex/build/matlab/Makefile.am

@@ -2,7 +2,7 @@ ACLOCAL_AMFLAGS = -I ../../../m4
 
 # libdynare++ must come before gensylv, k_order_perturbation, dynare_simul_
 if DO_SOMETHING
-SUBDIRS = mjdgges kronecker bytecode libdynare++ gensylv k_order_perturbation dynare_simul_
+SUBDIRS = mjdgges kronecker bytecode libdynare++ gensylv k_order_perturbation dynare_simul_ loglikelihood
 if HAVE_GSL
 SUBDIRS += swz
 endif

mex/build/matlab/configure.ac

@@ -102,6 +102,7 @@ AC_CONFIG_FILES([Makefile
                  bytecode/Makefile
                  k_order_perturbation/Makefile
                  dynare_simul_/Makefile
-                 swz/Makefile])
+                 swz/Makefile
+                 loglikelihood/Makefile])
 
 AC_OUTPUT

mex/build/matlab/loglikelihood/Makefile.am

+include ../mex.am
+include ../../loglikelihood.am

mex/build/octave/Makefile.am

@@ -2,7 +2,7 @@ ACLOCAL_AMFLAGS = -I ../../../m4
 
 # libdynare++ must come before gensylv, k_order_perturbation, dynare_simul_
 if DO_SOMETHING
-SUBDIRS = mjdgges kronecker bytecode libdynare++ gensylv k_order_perturbation dynare_simul_
+SUBDIRS = mjdgges kronecker bytecode libdynare++ gensylv k_order_perturbation dynare_simul_ loglikelihood
 if HAVE_GSL
 SUBDIRS += swz
 endif

mex/build/octave/configure.ac

@@ -84,6 +84,7 @@ AC_CONFIG_FILES([Makefile
                  gensylv/Makefile
                  k_order_perturbation/Makefile
                  dynare_simul_/Makefile
-                 swz/Makefile])
+                 swz/Makefile
+                 loglikelihood/Makefile])
 
 AC_OUTPUT

mex/build/octave/loglikelihood/Makefile.am

+include ../mex.am
+include ../../loglikelihood.am

mex/sources/estimation/Makefile.am

@@ -11,19 +11,25 @@ EXTRA_DIST = \
 	DecisionRules.hh \
 	DetrendData.cc \
 	DetrendData.hh \
+	EstimatedParameter.cc \
 	EstimatedParameter.hh \
+	EstimatedParametersDescription.cc \
 	EstimatedParametersDescription.hh \
+	EstimationSubsample.cc \
 	EstimationSubsample.hh \
 	InitializeKalmanFilter.cc \
 	InitializeKalmanFilter.hh \
 	KalmanFilter.cc \
 	KalmanFilter.hh \
-	LogLikelihoodMain.hh \
 	LogLikelihoodSubSample.cc \
 	LogLikelihoodSubSample.hh \
+	LogLikelihoodMain.hh \
+	LogLikelihoodMain.cc \
 	ModelSolution.cc \
 	ModelSolution.hh \
+	Prior.cc \
 	Prior.hh \
+	loglikelihood.cc \
 	utils/dynamic_dll.cc \
 	utils/dynamic_dll.hh \
 	utils/ts_exception.h

mex/sources/estimation/Prior.cc

@@ -36,3 +36,22 @@ Prior::~Prior()
 
 }
 
+Prior *
+Prior::constructPrior(pShape shape, double mean, double standard, double lower_bound, double upper_bound, double fhp, double shp)
+{
+  switch (shape)
+    {
+    case Beta:
+      return new BetaPrior(mean, standard, lower_bound, upper_bound, fhp, shp);
+    case Gamma:
+      return new GammaPrior(mean, standard, lower_bound, upper_bound, fhp, shp);
+    case Gaussian:
+      return new GaussianPrior(mean, standard, lower_bound, upper_bound, fhp, shp);
+    case Inv_gamma_1:
+      return new InvGamma1_Prior(mean, standard, lower_bound, upper_bound, fhp, shp);
+    case Uniform:
+      return new UniformPrior(mean, standard, lower_bound, upper_bound, fhp, shp);
+    case Inv_gamma_2:
+      return new InvGamma2_Prior(mean, standard, lower_bound, upper_bound, fhp, shp);
+    }
+}

mex/sources/estimation/Prior.hh

@@ -87,6 +87,7 @@ public:
     return 0.0;
   };
 
+  static Prior *constructPrior(pShape shape, double mean, double standard, double lower_bound, double upper_bound, double fhp, double shp);
 };
 
 struct BetaPrior : public Prior

mex/sources/estimation/libmat/Makefile.am

@@ -5,15 +5,18 @@ endif
 endif
 
 EXTRA_DIST = \
+	Matrix.hh \
+	Matrix.cc \
+	Vector.hh \
+	Vector.cc \
 	BlasBindings.hh \
 	DiscLyapFast.hh \
 	GeneralizedSchurDecomposition.cc \
 	GeneralizedSchurDecomposition.hh \
+	LapackBindings.hh \
 	LUSolver.cc \
 	LUSolver.hh \
-	Matrix.cc \
-	Matrix.hh \
 	QRDecomposition.cc \
 	QRDecomposition.hh \
-	Vector.cc \
-	Vector.hh
+	VDVEigDecomposition.cc \
+	VDVEigDecomposition.hh

mex/sources/estimation/loglikelihood.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 <string>
+#include <vector>
+#include <algorithm>
+#include <functional>
+
+#include "Vector.hh"
+#include "Matrix.hh"
+#include "LogLikelihoodMain.hh"
+
+#include "mex.h"
+
+void
+fillEstParamsInfo(const mxArray *estim_params_info, EstimatedParameter::pType type,
+                  std::vector<EstimatedParameter> &estParamsInfo)
+{
+  size_t m = mxGetM(estim_params_info), n = mxGetN(estim_params_info);
+  MatrixConstView epi(mxGetPr(estim_params_info), m, n, m);
+  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
+      double low_bound = epi(i, col++);
+      double up_bound = epi(i, col++);
+      Prior::pShape shape = (Prior::pShape) epi(i, col++);
+      double mean = epi(i, col++);
+      double std = epi(i, col++);
+      double p3 = epi(i, col++);
+      double p4 = epi(i, col++);
+
+      //      Prior *p = Prior::constructPrior(shape, mean, std, low_bound, up_bound, p3, p4);
+      Prior *p = NULL;
+
+      // Only one subsample
+      std::vector<size_t> subSampleIDs;
+      subSampleIDs.push_back(0);
+
+      estParamsInfo.push_back(EstimatedParameter(type, id1, id2, subSampleIDs,
+                                                 low_bound, up_bound, p));
+    }
+}
+
+double
+loglikelihood(const VectorConstView &estParams, const MatrixConstView &data,
+              const std::string &mexext)
+{
+  // Retrieve pointers to global variables
+  const mxArray *M_ = mexGetVariablePtr("global", "M_");
+  const mxArray *oo_ = mexGetVariablePtr("global", "oo_");
+  const mxArray *options_ = mexGetVariablePtr("global", "options_");
+  const mxArray *estim_params_ = mexGetVariablePtr("global", "estim_params_");
+
+  // Construct arguments of constructor of LogLikelihoodMain
+  char *fName = mxArrayToString(mxGetField(M_, 0, "fname"));
+  std::string dynamicDllFile(fName);
+  mxFree(fName);
+  dynamicDllFile += "_dynamic" + mexext;
+
+  size_t n_endo = (size_t) *mxGetPr(mxGetField(M_, 0, "endo_nbr"));
+  size_t n_exo = (size_t) *mxGetPr(mxGetField(M_, 0, "exo_nbr"));
+  size_t n_param = (size_t) *mxGetPr(mxGetField(M_, 0, "param_nbr"));
+  size_t n_estParams = estParams.getSize();
+
+  std::vector<size_t> zeta_fwrd, zeta_back, zeta_mixed, zeta_static;
+  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)
+    mexErrMsgTxt("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"));
+
+  std::vector<size_t> varobs;
+  const mxArray *varobs_mx = mxGetField(options_, 0, "varobs_id");
+  if (mxGetM(varobs_mx) != 1)
+    mexErrMsgTxt("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)
+    mexErrMsgTxt("Data has not 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);
+
+  // Allocate LogLikelihoodMain object
+  int info;
+  LogLikelihoodMain llm(dynamicDllFile, epd, n_endo, n_exo, zeta_fwrd, zeta_back, zeta_mixed, zeta_static,
+                        qz_criterium, varobs, riccati_tol, lyapunov_tol, info);
+
+  // Construct arguments of compute() method
+  Matrix steadyState(n_endo, 1);
+  mat::get_col(steadyState, 0) = VectorConstView(mxGetPr(mxGetField(oo_, 0, "steady_state")), n_endo, 1);
+
+  Vector estParams2(n_estParams);
+  estParams2 = estParams;
+  Vector deepParams(n_param);
+  deepParams = VectorConstView(mxGetPr(mxGetField(M_, 0, "params")), n_param, 1);
+  Matrix Q(n_exo);
+  Q = MatrixConstView(mxGetPr(mxGetField(M_, 0, "Sigma_e")), n_exo, n_exo, 1);
+  Matrix H(n_varobs);
+  const mxArray *H_mx = mxGetField(M_, 0, "H");
+  if (mxGetM(H_mx) == 1 && mxGetN(H_mx) == 1 && *mxGetPr(H_mx) == 0)
+    H.setAll(0.0);
+  else
+    H = MatrixConstView(mxGetPr(mxGetField(M_, 0, "H")), n_varobs, n_varobs, 1);
+
+  // Compute the likelihood
+  double lik = llm.compute(steadyState, estParams2, deepParams, data, Q, H, 0, info);
+
+  // Cleanups
+  /*
+  for (std::vector<EstimatedParameter>::iterator it = estParamsInfo.begin();
+       it != estParamsInfo.end(); it++)
+    delete it->prior;
+  */
+
+  return lik;
+}
+
+void
+mexFunction(int nlhs, mxArray *plhs[],
+            int nrhs, const mxArray *prhs[])
+{
+  if (nrhs != 3)
+    mexErrMsgTxt("loglikelihood: exactly three arguments are required.");
+  if (nlhs != 1)
+    mexErrMsgTxt("loglikelihood: exactly one return argument is required.");
+
+  // Check and retrieve the arguments
+
+  if (!mxIsDouble(prhs[0]) || mxGetN(prhs[0]) != 1)
+    mexErrMsgTxt("First argument must be a column vector of double-precision numbers");
+
+  VectorConstView estParams(mxGetPr(prhs[0]), mxGetM(prhs[0]), 1);
+
+  if (!mxIsDouble(prhs[1]))
+    mexErrMsgTxt("Second argument must be a matrix of double-precision numbers");
+
+  MatrixConstView data(mxGetPr(prhs[1]), mxGetM(prhs[1]), mxGetN(prhs[1]), mxGetM(prhs[1]));
+
+  if (!mxIsChar(prhs[2]))
+    mexErrMsgTxt("Third argument must be a character string");
+
+  char *mexext_mx = mxArrayToString(prhs[2]);
+  std::string mexext(mexext_mx);
+  mxFree(mexext_mx);
+
+  // Compute and return the value
+  double lik = loglikelihood(estParams, data, mexext);
+
+  plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL); 
+  *mxGetPr(plhs[0]) = lik;
+}

preprocessor/SymbolTable.cc

@@ -250,6 +250,12 @@ SymbolTable::writeOutput(ostream &output) const throw (NotYetFrozenException)
       for (vector<int>::const_iterator it = varobs.begin();
            it != varobs.end(); it++)
         output << "options_.varobs = strvcat(options_.varobs, '" << getName(*it) << "');" << endl;
+
+      output << "options_.varobs_id = [ ";
+      for (vector<int>::const_iterator it = varobs.begin();
+           it != varobs.end(); it++)
+        output << getTypeSpecificID(*it)+1 << " ";
+      output << " ];"  << endl;
     }
 }