From 75fbd3852434dcdd9db12b6a55cf659ae5d632a5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= Date: Thu, 7 Feb 2013 18:33:55 +0100 Subject: [PATCH] Handle constant in Kalman filter Previously, the filter was only working with options_.noconstant, i.e. when the data are centered. --- mex/sources/estimation/DetrendData.cc | 13 +++++++++++-- mex/sources/estimation/DetrendData.hh | 6 +++++- mex/sources/estimation/InitializeKalmanFilter.cc | 6 +++++- mex/sources/estimation/InitializeKalmanFilter.hh | 4 +++- mex/sources/estimation/KalmanFilter.cc | 5 +++-- mex/sources/estimation/KalmanFilter.hh | 3 ++- mex/sources/estimation/LogLikelihoodMain.cc | 7 ++++--- mex/sources/estimation/LogLikelihoodMain.hh | 5 +++-- mex/sources/estimation/LogLikelihoodSubSample.cc | 4 ++-- mex/sources/estimation/LogLikelihoodSubSample.hh | 2 +- mex/sources/estimation/LogPosteriorDensity.cc | 7 ++++--- mex/sources/estimation/LogPosteriorDensity.hh | 5 +++-- mex/sources/estimation/logMHMCMCposterior.cc | 5 ++++- mex/sources/estimation/logposterior.cc | 4 +++- 14 files changed, 53 insertions(+), 23 deletions(-) diff --git a/mex/sources/estimation/DetrendData.cc b/mex/sources/estimation/DetrendData.cc index 5d2a63506..6349932c9 100644 --- a/mex/sources/estimation/DetrendData.cc +++ b/mex/sources/estimation/DetrendData.cc @@ -25,7 +25,9 @@ #include "DetrendData.hh" -DetrendData::DetrendData() +DetrendData::DetrendData(const std::vector &varobs_arg, + bool noconstant_arg) + : varobs(varobs_arg), noconstant(noconstant_arg) { }; @@ -33,6 +35,13 @@ void DetrendData::detrend(const VectorView &SteadyState, const MatrixConstView &dataView, MatrixView &detrendedDataView) { - detrendedDataView = dataView; + if (noconstant) + detrendedDataView = dataView; + else + { + for (size_t i = 0; i < varobs.size(); i++) + for (size_t j = 0; j < dataView.getCols(); j++) + detrendedDataView(i, j) = dataView(i, j) - SteadyState(varobs[i]); + } }; diff --git a/mex/sources/estimation/DetrendData.hh b/mex/sources/estimation/DetrendData.hh index 1e91abc6f..59cfc3d0f 100644 --- a/mex/sources/estimation/DetrendData.hh +++ b/mex/sources/estimation/DetrendData.hh @@ -32,9 +32,13 @@ class DetrendData public: virtual ~DetrendData(){}; - DetrendData(); + DetrendData(const std::vector &varobs_arg, bool noconstant_arg); void detrend(const VectorView &SteadyState, const MatrixConstView &dataView, MatrixView &detrendedDataView); +private: + const std::vector varobs; + const bool noconstant; + }; #endif // !defined(312823A1_6248_4af0_B204_DB22F1237E9B__INCLUDED_) diff --git a/mex/sources/estimation/InitializeKalmanFilter.cc b/mex/sources/estimation/InitializeKalmanFilter.cc index 7ed474996..ec9b93485 100644 --- a/mex/sources/estimation/InitializeKalmanFilter.cc +++ b/mex/sources/estimation/InitializeKalmanFilter.cc @@ -34,10 +34,14 @@ InitializeKalmanFilter::InitializeKalmanFilter(const std::string &dynamicDllFile const std::vector &zeta_fwrd_arg, const std::vector &zeta_back_arg, const std::vector &zeta_mixed_arg, const std::vector &zeta_static_arg, const std::vector &zeta_varobs_back_mixed_arg, + const std::vector &varobs_arg, double qz_criterium_arg, - double lyapunov_tol_arg, int &info) : + double lyapunov_tol_arg, + bool noconstant_arg, + int &info) : lyapunov_tol(lyapunov_tol_arg), zeta_varobs_back_mixed(zeta_varobs_back_mixed_arg), + detrendData(varobs_arg, noconstant_arg), modelSolution(dynamicDllFile, n_endo_arg, n_exo_arg, zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg, zeta_static_arg, qz_criterium_arg), discLyapFast(zeta_varobs_back_mixed.size()), diff --git a/mex/sources/estimation/InitializeKalmanFilter.hh b/mex/sources/estimation/InitializeKalmanFilter.hh index 4ddced4de..0a4414794 100644 --- a/mex/sources/estimation/InitializeKalmanFilter.hh +++ b/mex/sources/estimation/InitializeKalmanFilter.hh @@ -46,7 +46,9 @@ public: InitializeKalmanFilter(const std::string &dynamicDllFile, size_t n_endo, size_t n_exo, const std::vector &zeta_fwrd_arg, const std::vector &zeta_back_arg, const std::vector &zeta_mixed_arg, const std::vector &zeta_static_arg, const std::vector &zeta_varobs_back_mixed_arg, - double qz_criterium_arg, double lyapunov_tol_arg, int &info); + const std::vector &varobs_arg, + double qz_criterium_arg, double lyapunov_tol_arg, + bool noconstant_arg, int &info); virtual ~InitializeKalmanFilter(); // initialise parameter dependent KF matrices only but not Ps template diff --git a/mex/sources/estimation/KalmanFilter.cc b/mex/sources/estimation/KalmanFilter.cc index 5b416e1e9..313f6b4ce 100644 --- a/mex/sources/estimation/KalmanFilter.cc +++ b/mex/sources/estimation/KalmanFilter.cc @@ -35,7 +35,8 @@ KalmanFilter::KalmanFilter(const std::string &dynamicDllFile, size_t n_endo, siz const std::vector &zeta_fwrd_arg, const std::vector &zeta_back_arg, const std::vector &zeta_mixed_arg, const std::vector &zeta_static_arg, double qz_criterium_arg, const std::vector &varobs_arg, - double riccati_tol_arg, double lyapunov_tol_arg, int &info) : + double riccati_tol_arg, double lyapunov_tol_arg, + bool noconstant_arg, int &info) : zeta_varobs_back_mixed(compute_zeta_varobs_back_mixed(zeta_back_arg, zeta_mixed_arg, varobs_arg)), Z(varobs_arg.size(), zeta_varobs_back_mixed.size()), Zt(Z.getCols(), Z.getRows()), T(zeta_varobs_back_mixed.size()), R(zeta_varobs_back_mixed.size(), n_exo), Pstar(zeta_varobs_back_mixed.size(), zeta_varobs_back_mixed.size()), Pinf(zeta_varobs_back_mixed.size(), zeta_varobs_back_mixed.size()), @@ -44,7 +45,7 @@ KalmanFilter::KalmanFilter(const std::string &dynamicDllFile, size_t n_endo, siz oldKFinv(zeta_varobs_back_mixed.size(), varobs_arg.size()), a_init(zeta_varobs_back_mixed.size()), a_new(zeta_varobs_back_mixed.size()), vt(varobs_arg.size()), vtFinv(varobs_arg.size()), riccati_tol(riccati_tol_arg), initKalmanFilter(dynamicDllFile, n_endo, n_exo, zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg, - zeta_static_arg, zeta_varobs_back_mixed, qz_criterium_arg, lyapunov_tol_arg, info), + zeta_static_arg, zeta_varobs_back_mixed, varobs_arg, qz_criterium_arg, lyapunov_tol_arg, noconstant_arg, info), FUTP(varobs_arg.size()*(varobs_arg.size()+1)/2) { Z.setAll(0.0); diff --git a/mex/sources/estimation/KalmanFilter.hh b/mex/sources/estimation/KalmanFilter.hh index be3c2e868..0a3ee50eb 100644 --- a/mex/sources/estimation/KalmanFilter.hh +++ b/mex/sources/estimation/KalmanFilter.hh @@ -52,7 +52,8 @@ public: KalmanFilter(const std::string &dynamicDllFile, size_t n_endo, size_t n_exo, const std::vector &zeta_fwrd_arg, const std::vector &zeta_back_arg, const std::vector &zeta_mixed_arg, const std::vector &zeta_static_arg, double qz_criterium_arg, const std::vector &varobs_arg, - double riccati_tol_arg, double lyapunov_tol_arg, int &info); + double riccati_tol_arg, double lyapunov_tol_arg, + bool noconstant_arg, int &info); template double compute(const MatrixConstView &dataView, Vec1 &steadyState, diff --git a/mex/sources/estimation/LogLikelihoodMain.cc b/mex/sources/estimation/LogLikelihoodMain.cc index 47318726d..a30de15d0 100644 --- a/mex/sources/estimation/LogLikelihoodMain.cc +++ b/mex/sources/estimation/LogLikelihoodMain.cc @@ -1,5 +1,5 @@ /* - * Copyright (C) 2009-2012 Dynare Team + * Copyright (C) 2009-2013 Dynare Team * * This file is part of Dynare. * @@ -28,11 +28,12 @@ LogLikelihoodMain::LogLikelihoodMain(const std::string &dynamicDllFile, EstimatedParametersDescription &estiParDesc, size_t n_endo, size_t n_exo, const std::vector &zeta_fwrd_arg, const std::vector &zeta_back_arg, const std::vector &zeta_mixed_arg, const std::vector &zeta_static_arg, const double qz_criterium, - const std::vector &varobs, double riccati_tol, double lyapunov_tol, int &info_arg) + const std::vector &varobs, double riccati_tol, double lyapunov_tol, + bool noconstant_arg, int &info_arg) : estSubsamples(estiParDesc.estSubsamples), logLikelihoodSubSample(dynamicDllFile, estiParDesc, n_endo, n_exo, zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg, zeta_static_arg, qz_criterium, - varobs, riccati_tol, lyapunov_tol, info_arg), + varobs, riccati_tol, lyapunov_tol, noconstant_arg, info_arg), vll(estiParDesc.getNumberOfPeriods()), // time dimension size of data detrendedData(varobs.size(), estiParDesc.getNumberOfPeriods()) { diff --git a/mex/sources/estimation/LogLikelihoodMain.hh b/mex/sources/estimation/LogLikelihoodMain.hh index 88a2b5bfa..3811104b2 100644 --- a/mex/sources/estimation/LogLikelihoodMain.hh +++ b/mex/sources/estimation/LogLikelihoodMain.hh @@ -1,5 +1,5 @@ /* - * Copyright (C) 2009-2012 Dynare Team + * Copyright (C) 2009-2013 Dynare Team * * This file is part of Dynare. * @@ -40,7 +40,8 @@ public: LogLikelihoodMain(const std::string &dynamicDllFile, EstimatedParametersDescription &estiParDesc, size_t n_endo, size_t n_exo, const std::vector &zeta_fwrd_arg, const std::vector &zeta_back_arg, const std::vector &zeta_mixed_arg, const std::vector &zeta_static_arg, const double qz_criterium_arg, const std::vector &varobs_arg, - double riccati_tol_arg, double lyapunov_tol_arg, int &info); + double riccati_tol_arg, double lyapunov_tol_arg, + bool noconstant_arg, int &info); /** * Compute method Inputs: diff --git a/mex/sources/estimation/LogLikelihoodSubSample.cc b/mex/sources/estimation/LogLikelihoodSubSample.cc index b84fa0d80..e16b8ebe9 100644 --- a/mex/sources/estimation/LogLikelihoodSubSample.cc +++ b/mex/sources/estimation/LogLikelihoodSubSample.cc @@ -32,10 +32,10 @@ LogLikelihoodSubSample::~LogLikelihoodSubSample() LogLikelihoodSubSample::LogLikelihoodSubSample(const std::string &dynamicDllFile, EstimatedParametersDescription &INestiParDesc, size_t n_endo, size_t n_exo, const std::vector &zeta_fwrd_arg, const std::vector &zeta_back_arg, const std::vector &zeta_mixed_arg, const std::vector &zeta_static_arg, const double qz_criterium, - const std::vector &varobs, double riccati_tol, double lyapunov_tol, int &INinfo) : + const std::vector &varobs, double riccati_tol, double lyapunov_tol, bool noconstant_arg, int &INinfo) : startPenalty(-1e8), estiParDesc(INestiParDesc), kalmanFilter(dynamicDllFile, n_endo, n_exo, zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg, zeta_static_arg, qz_criterium, - varobs, riccati_tol, lyapunov_tol, INinfo), eigQ(n_exo), eigH(varobs.size()), info(INinfo) + varobs, riccati_tol, lyapunov_tol, noconstant_arg, INinfo), eigQ(n_exo), eigH(varobs.size()), info(INinfo) { }; diff --git a/mex/sources/estimation/LogLikelihoodSubSample.hh b/mex/sources/estimation/LogLikelihoodSubSample.hh index 80f15489b..8be88e0fc 100644 --- a/mex/sources/estimation/LogLikelihoodSubSample.hh +++ b/mex/sources/estimation/LogLikelihoodSubSample.hh @@ -39,7 +39,7 @@ public: LogLikelihoodSubSample(const std::string &dynamicDllFile, EstimatedParametersDescription &estiParDesc, size_t n_endo, size_t n_exo, const std::vector &zeta_fwrd_arg, const std::vector &zeta_back_arg, const std::vector &zeta_mixed_arg, const std::vector &zeta_static_arg, const double qz_criterium, - const std::vector &varobs_arg, double riccati_tol_in, double lyapunov_tol, int &info); + const std::vector &varobs_arg, double riccati_tol_in, double lyapunov_tol, bool noconstant_arg, int &info); template double compute(VEC1 &steadyState, const MatrixConstView &dataView, VEC2 &estParams, VectorView &deepParams, diff --git a/mex/sources/estimation/LogPosteriorDensity.cc b/mex/sources/estimation/LogPosteriorDensity.cc index d938b9565..114be39ea 100644 --- a/mex/sources/estimation/LogPosteriorDensity.cc +++ b/mex/sources/estimation/LogPosteriorDensity.cc @@ -1,5 +1,5 @@ /* - * Copyright (C) 2009-2012 Dynare Team + * Copyright (C) 2009-2013 Dynare Team * * This file is part of Dynare. * @@ -32,10 +32,11 @@ LogPosteriorDensity::~LogPosteriorDensity() LogPosteriorDensity::LogPosteriorDensity(const std::string &modName, EstimatedParametersDescription &estParamsDesc, size_t n_endo, size_t n_exo, const std::vector &zeta_fwrd_arg, const std::vector &zeta_back_arg, const std::vector &zeta_mixed_arg, const std::vector &zeta_static_arg, const double qz_criterium_arg, const std::vector &varobs_arg, - double riccati_tol_arg, double lyapunov_tol_arg, int &info_arg) : + double riccati_tol_arg, double lyapunov_tol_arg, + bool noconstant_arg, int &info_arg) : logPriorDensity(estParamsDesc), logLikelihoodMain(modName, estParamsDesc, n_endo, n_exo, zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg, - zeta_static_arg, qz_criterium_arg, varobs_arg, riccati_tol_arg, lyapunov_tol_arg, info_arg) + zeta_static_arg, qz_criterium_arg, varobs_arg, riccati_tol_arg, lyapunov_tol_arg, noconstant_arg, info_arg) { } diff --git a/mex/sources/estimation/LogPosteriorDensity.hh b/mex/sources/estimation/LogPosteriorDensity.hh index 4ff4a0d59..e11795c9e 100644 --- a/mex/sources/estimation/LogPosteriorDensity.hh +++ b/mex/sources/estimation/LogPosteriorDensity.hh @@ -1,5 +1,5 @@ /* - * Copyright (C) 2009-2012 Dynare Team + * Copyright (C) 2009-2013 Dynare Team * * This file is part of Dynare. * @@ -47,7 +47,8 @@ public: LogPosteriorDensity(const std::string &modName, EstimatedParametersDescription &estParamsDesc, size_t n_endo, size_t n_exo, const std::vector &zeta_fwrd_arg, const std::vector &zeta_back_arg, const std::vector &zeta_mixed_arg, const std::vector &zeta_static_arg, const double qz_criterium_arg, const std::vector &varobs_arg, - double riccati_tol_arg, double lyapunov_tol_arg, int &info_arg); + double riccati_tol_arg, double lyapunov_tol_arg, + bool noconstant_arg, int &info_arg); template double diff --git a/mex/sources/estimation/logMHMCMCposterior.cc b/mex/sources/estimation/logMHMCMCposterior.cc index cd97f621a..7cdef7e31 100644 --- a/mex/sources/estimation/logMHMCMCposterior.cc +++ b/mex/sources/estimation/logMHMCMCposterior.cc @@ -711,10 +711,13 @@ logMCMCposterior(VectorConstView &estParams, const MatrixConstView &data, estParamsInfo); EstimatedParametersDescription epd(estSubsamples, estParamsInfo); + + bool noconstant = (bool) *mxGetPr(mxGetField(options_, 0, "noconstant")); + // Allocate LogPosteriorDensity object int info; LogPosteriorDensity lpd(dynamicDllFile, epd, n_endo, n_exo, zeta_fwrd, zeta_back, zeta_mixed, zeta_static, - qz_criterium, varobs, riccati_tol, lyapunov_tol, info); + qz_criterium, varobs, riccati_tol, lyapunov_tol, noconstant, info); // Construct MHMCMC Sampler RandomWalkMetropolisHastings rwmh(estParams.getSize()); diff --git a/mex/sources/estimation/logposterior.cc b/mex/sources/estimation/logposterior.cc index f0910ca3e..ca370dca0 100644 --- a/mex/sources/estimation/logposterior.cc +++ b/mex/sources/estimation/logposterior.cc @@ -181,9 +181,11 @@ logposterior(VEC1 &estParams, const MatrixConstView &data, EstimatedParametersDescription epd(estSubsamples, estParamsInfo); + bool noconstant = (bool) *mxGetPr(mxGetField(options_, 0, "noconstant")); + // Allocate LogPosteriorDensity object LogPosteriorDensity lpd(dynamicDllFile, epd, n_endo, n_exo, zeta_fwrd, zeta_back, zeta_mixed, zeta_static, - qz_criterium, varobs, riccati_tol, lyapunov_tol, info); + qz_criterium, varobs, riccati_tol, lyapunov_tol, noconstant, info); // Construct arguments of compute() method