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

@@ -25,7 +25,9 @@
 
 #include "DetrendData.hh"
 
-DetrendData::DetrendData()
+DetrendData::DetrendData(const std::vector<size_t> &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]);
+    }
 };
 

mex/sources/estimation/DetrendData.hh

@@ -32,9 +32,13 @@ class DetrendData
 
 public:
   virtual ~DetrendData(){};
-  DetrendData();
+  DetrendData(const std::vector<size_t> &varobs_arg, bool noconstant_arg);
   void detrend(const VectorView &SteadyState, const MatrixConstView &dataView, MatrixView &detrendedDataView);
 
+private:
+  const std::vector<size_t> varobs;
+  const bool noconstant;
+
 };
 
 #endif // !defined(312823A1_6248_4af0_B204_DB22F1237E9B__INCLUDED_)

mex/sources/estimation/InitializeKalmanFilter.cc

@@ -34,10 +34,14 @@ InitializeKalmanFilter::InitializeKalmanFilter(const std::string &dynamicDllFile
                                                const std::vector<size_t> &zeta_fwrd_arg, const std::vector<size_t> &zeta_back_arg,
                                                const std::vector<size_t> &zeta_mixed_arg, const std::vector<size_t> &zeta_static_arg,
                                                const std::vector<size_t> &zeta_varobs_back_mixed_arg,
+                                               const std::vector<size_t> &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()),

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<size_t> &zeta_fwrd_arg,
                          const std::vector<size_t> &zeta_back_arg, const std::vector<size_t> &zeta_mixed_arg, const std::vector<size_t> &zeta_static_arg,
                          const std::vector<size_t> &zeta_varobs_back_mixed_arg,
-                         double qz_criterium_arg, double lyapunov_tol_arg, int &info);
+                         const std::vector<size_t> &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 <class Vec1, class Vec2, class Mat1, class Mat2>

mex/sources/estimation/KalmanFilter.cc

@@ -35,7 +35,8 @@ KalmanFilter::KalmanFilter(const std::string &dynamicDllFile, size_t n_endo, siz
                            const std::vector<size_t> &zeta_fwrd_arg, const std::vector<size_t> &zeta_back_arg,
                            const std::vector<size_t> &zeta_mixed_arg, const std::vector<size_t> &zeta_static_arg,
                            double qz_criterium_arg, const std::vector<size_t> &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);

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<size_t> &zeta_fwrd_arg,
                const std::vector<size_t> &zeta_back_arg, const std::vector<size_t> &zeta_mixed_arg, const std::vector<size_t> &zeta_static_arg,
                double qz_criterium_arg, const std::vector<size_t> &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 <class Vec1, class Vec2, class Mat1>
   double compute(const MatrixConstView &dataView, Vec1 &steadyState,

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<size_t> &zeta_fwrd_arg, const std::vector<size_t> &zeta_back_arg,
                                      const std::vector<size_t> &zeta_mixed_arg, const std::vector<size_t> &zeta_static_arg, const double qz_criterium,
-                                     const std::vector<size_t> &varobs, double riccati_tol, double lyapunov_tol, int &info_arg)
+                                     const std::vector<size_t> &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())
 {

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<size_t> &zeta_fwrd_arg, const std::vector<size_t> &zeta_back_arg, const std::vector<size_t> &zeta_mixed_arg,
                     const std::vector<size_t> &zeta_static_arg, const double qz_criterium_arg, const std::vector<size_t> &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:

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<size_t> &zeta_fwrd_arg, const std::vector<size_t> &zeta_back_arg,
                                                const std::vector<size_t> &zeta_mixed_arg, const std::vector<size_t> &zeta_static_arg, const double qz_criterium,
-                                               const std::vector<size_t> &varobs, double riccati_tol, double lyapunov_tol, int &INinfo) :
+                                               const std::vector<size_t> &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)
 {
 };
 

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<size_t> &zeta_fwrd_arg, const std::vector<size_t> &zeta_back_arg,
                          const std::vector<size_t> &zeta_mixed_arg, const std::vector<size_t> &zeta_static_arg, const double qz_criterium,
-                         const std::vector<size_t> &varobs_arg, double riccati_tol_in, double lyapunov_tol, int &info);
+                         const std::vector<size_t> &varobs_arg, double riccati_tol_in, double lyapunov_tol, bool noconstant_arg, int &info);
 
   template <class VEC1, class VEC2>
   double compute(VEC1 &steadyState, const MatrixConstView &dataView, VEC2 &estParams, VectorView &deepParams,

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<size_t> &zeta_fwrd_arg, const std::vector<size_t> &zeta_back_arg, const std::vector<size_t> &zeta_mixed_arg,
                                          const std::vector<size_t> &zeta_static_arg, const double qz_criterium_arg, const std::vector<size_t> &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)
 {
 
 }

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<size_t> &zeta_fwrd_arg, const std::vector<size_t> &zeta_back_arg, const std::vector<size_t> &zeta_mixed_arg,
                       const std::vector<size_t> &zeta_static_arg, const double qz_criterium_arg, const std::vector<size_t> &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 <class VEC1, class VEC2>
   double

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());

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