diff --git a/mex/build/loglikelihood.am b/mex/build/loglikelihood.am new file mode 100644 index 000000000..4664b4bf2 --- /dev/null +++ b/mex/build/loglikelihood.am @@ -0,0 +1,53 @@ +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 diff --git a/mex/build/matlab/Makefile.am b/mex/build/matlab/Makefile.am index a82facefc..f823b41ea 100644 --- a/mex/build/matlab/Makefile.am +++ b/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 diff --git a/mex/build/matlab/configure.ac b/mex/build/matlab/configure.ac index f1e77acaa..e39f90c56 100644 --- a/mex/build/matlab/configure.ac +++ b/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 diff --git a/mex/build/matlab/loglikelihood/Makefile.am b/mex/build/matlab/loglikelihood/Makefile.am new file mode 100644 index 000000000..c59b3ffd4 --- /dev/null +++ b/mex/build/matlab/loglikelihood/Makefile.am @@ -0,0 +1,2 @@ +include ../mex.am +include ../../loglikelihood.am diff --git a/mex/build/octave/Makefile.am b/mex/build/octave/Makefile.am index 4a3c8e66c..4bf868a4f 100644 --- a/mex/build/octave/Makefile.am +++ b/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 diff --git a/mex/build/octave/configure.ac b/mex/build/octave/configure.ac index 580ab6965..626ea6eda 100644 --- a/mex/build/octave/configure.ac +++ b/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 diff --git a/mex/build/octave/loglikelihood/Makefile.am b/mex/build/octave/loglikelihood/Makefile.am new file mode 100644 index 000000000..c59b3ffd4 --- /dev/null +++ b/mex/build/octave/loglikelihood/Makefile.am @@ -0,0 +1,2 @@ +include ../mex.am +include ../../loglikelihood.am diff --git a/mex/sources/estimation/Makefile.am b/mex/sources/estimation/Makefile.am index e820c961e..930895f2b 100644 --- a/mex/sources/estimation/Makefile.am +++ b/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 diff --git a/mex/sources/estimation/Prior.cc b/mex/sources/estimation/Prior.cc index 1cc47dec5..1272a39c0 100644 --- a/mex/sources/estimation/Prior.cc +++ b/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); + } +} diff --git a/mex/sources/estimation/Prior.hh b/mex/sources/estimation/Prior.hh index 37bc72ab0..64a7e6161 100644 --- a/mex/sources/estimation/Prior.hh +++ b/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 diff --git a/mex/sources/estimation/libmat/Makefile.am b/mex/sources/estimation/libmat/Makefile.am index ef4a5a450..b73cf569f 100644 --- a/mex/sources/estimation/libmat/Makefile.am +++ b/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 diff --git a/mex/sources/estimation/loglikelihood.cc b/mex/sources/estimation/loglikelihood.cc new file mode 100644 index 000000000..93b0e84ff --- /dev/null +++ b/mex/sources/estimation/loglikelihood.cc @@ -0,0 +1,204 @@ +/* + * 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 . + */ + +#include +#include +#include +#include + +#include "Vector.hh" +#include "Matrix.hh" +#include "LogLikelihoodMain.hh" + +#include "mex.h" + +void +fillEstParamsInfo(const mxArray *estim_params_info, EstimatedParameter::pType type, + std::vector &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 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 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 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(), 1)); + + if (data.getRows() != n_varobs) + mexErrMsgTxt("Data has not as many rows as there are observed variables"); + + std::vector estSubsamples; + estSubsamples.push_back(EstimationSubsample(0, data.getCols() - 1)); + + std::vector 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::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; +} diff --git a/preprocessor/SymbolTable.cc b/preprocessor/SymbolTable.cc index ae4af72c3..b3e4ad12f 100644 --- a/preprocessor/SymbolTable.cc +++ b/preprocessor/SymbolTable.cc @@ -250,6 +250,12 @@ SymbolTable::writeOutput(ostream &output) const throw (NotYetFrozenException) for (vector::const_iterator it = varobs.begin(); it != varobs.end(); it++) output << "options_.varobs = strvcat(options_.varobs, '" << getName(*it) << "');" << endl; + + output << "options_.varobs_id = [ "; + for (vector::const_iterator it = varobs.begin(); + it != varobs.end(); it++) + output << getTypeSpecificID(*it)+1 << " "; + output << " ];" << endl; } }