Estimation C++DLL: First draft for LogPosteriorDensity, LogPriorDensity and...

Estimation C++DLL: First draft for LogPosteriorDensity, LogPriorDensity and associated changes: random number generation still missing

mex/sources/estimation/EstimatedParameter.cc

@@ -27,7 +27,7 @@
 
 EstimatedParameter::EstimatedParameter(const EstimatedParameter::pType type_arg,
                                        size_t ID1_arg, size_t ID2_arg, const std::vector<size_t> &subSampleIDs_arg,
-                                       double lower_bound_arg, double upper_bound_arg, Prior prior_arg) :
+                                       double lower_bound_arg, double upper_bound_arg, Prior* prior_arg) :
   ptype(type_arg), ID1(ID1_arg), ID2(ID2_arg),
   lower_bound(lower_bound_arg), upper_bound(upper_bound_arg), prior(prior_arg),
   subSampleIDs(subSampleIDs_arg)

mex/sources/estimation/EstimatedParameter.hh

@@ -43,7 +43,7 @@ public:
 
   EstimatedParameter(const EstimatedParameter::pType type,
                      size_t ID1, size_t ID2, const std::vector<size_t> &subSampleIDs,
-                     double lower_bound, double upper_bound, Prior prior
+                     double lower_bound, double upper_bound, Prior* prior
                      );
   virtual ~EstimatedParameter();
 
@@ -52,7 +52,7 @@ public:
   const size_t ID2;
   const double lower_bound;
   const double upper_bound;
-  const Prior prior;
+  Prior* prior;
   const std::vector<size_t> subSampleIDs;
 };
 

mex/sources/estimation/LogLikelihoodMain.hh

@@ -50,7 +50,8 @@ public:
    * Q and H KF matrices of shock and measurement error varinaces and covariances
    * KF logLikelihood calculation start period.
    */
-  double compute(Matrix &steadyState, const Vector &estParams, Vector &deepParams, const MatrixConstView &data, Matrix &Q, Matrix &H, size_t presampleStart, int &info); // for calls from estimation and to set Steady State
+  double compute(Matrix &steadyState, const Vector &estParams, Vector &deepParams, const MatrixConstView &data, 
+    Matrix &Q, Matrix &H, size_t presampleStart, int &info); // for calls from estimation and to set Steady State
 
   Vector &
   getVll()

mex/sources/estimation/LogLikelihoodSubSample.hh

@@ -26,8 +26,8 @@
 #if !defined(DF8B7AF5_8169_4587_9037_2CD2C82E2DDF__INCLUDED_)
 #define DF8B7AF5_8169_4587_9037_2CD2C82E2DDF__INCLUDED_
 
-#include "KalmanFilter.hh"
 #include "EstimatedParametersDescription.hh"
+#include "KalmanFilter.hh"
 #include "VDVEigDecomposition.hh"
 
 class LogLikelihoodSubSample {

mex/sources/estimation/LogPosteriorDensity.cc

+/*
+ * Copyright (C) 2009-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/>.
+ */
+
+///////////////////////////////////////////////////////////
+//  LogPosteriorDensity.cpp
+//  Implementation of the Class LogPosteriorDensity
+//  Created on:      10-Feb-2010 20:54:18
+///////////////////////////////////////////////////////////
+
+#include "LogPosteriorDensity.hh"
+
+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) :
+  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)
+{
+
+}
+
+double
+LogPosteriorDensity::compute(Matrix &steadyState, const Vector &estParams, Vector &deepParams, const MatrixConstView &data, Matrix &Q, Matrix &H, size_t presampleStart, int &info)
+{
+  return logLikelihoodMain.compute(steadyState, estParams, deepParams, data, Q, H, presampleStart, info)
+         +logPriorDensity.compute(estParams);
+}
+
+/**
+ * vector of log likelihoods for each Kalman step
+ */
+Vector &
+LogPosteriorDensity::getLikVector()
+{
+  return logLikelihoodMain.getVll();
+}
+
+/**
+ * log likelihood as summ of Vll for each Kalman step
+ */
+double
+LogPosteriorDensity::getLogPosteriorDensity()
+{
+  return logLikelihoodMain.getLogLikelihood()+logPriorDensity.getLogPriorDensity();
+}
+

mex/sources/estimation/LogPosteriorDensity.hh

+/*
+ * Copyright (C) 2009-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/>.
+ */
+
+///////////////////////////////////////////////////////////
+//  LogPosteriorDensity.hh
+//  Implementation of the Class LogPosteriorDensity
+//  Created on:      10-Feb-2010 20:54:18
+///////////////////////////////////////////////////////////
+
+#if !defined(LPD_052A31B5_53BF_4904_AD80_863B52827973__INCLUDED_)
+#define LPD_052A31B5_53BF_4904_AD80_863B52827973__INCLUDED_
+
+#include "EstimatedParametersDescription.hh"
+#include "LogPriorDensity.hh"
+#include "LogLikelihoodMain.hh"
+
+/**
+ * Class that calculates Log Posterior Density using kalman, based on Dynare
+ * DsgeLikelihood.m
+ */
+class LogPosteriorDensity {
+
+private:
+  LogPriorDensity logPriorDensity;
+  LogLikelihoodMain logLikelihoodMain;
+  double logPosteriorDensity;
+
+public:
+  virtual ~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 compute(Matrix &steadyState, const Vector &estParams, Vector &deepParams,
+                 const MatrixConstView &data, Matrix &Q, Matrix &H, size_t presampleStart, int &info);
+  Vector&getLikVector();
+  double getLogPosteriorDensity();
+
+};
+
+#endif // !defined(052A31B5_53BF_4904_AD80_863B52827973__INCLUDED_)

mex/sources/estimation/LogPriorDensity.cc

+/*
+ * Copyright (C) 2009-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/>.
+ */
+
+///////////////////////////////////////////////////////////
+//  LogPriorDensity.cpp
+//  Implementation of the Class LogPriorDensity
+//  Created on:      10-Feb-2010 20:56:08
+///////////////////////////////////////////////////////////
+
+#include "LogPriorDensity.hh"
+
+LogPriorDensity::~LogPriorDensity()
+{
+};
+
+LogPriorDensity::LogPriorDensity(EstimatedParametersDescription &estParsDesc_arg) :
+  estParsDesc(estParsDesc_arg)
+{
+};
+
+double
+LogPriorDensity::compute(const Vector &ep)
+{
+  assert(estParsDesc.estParams.size() == ep.getSize());
+
+  for (size_t i = 0; i <  ep.getSize(); ++i)
+    {
+      logPriorDensity += ((*(estParsDesc.estParams[i]).prior)).pdf(ep(i));
+      if (std::isinf(abs(logPriorDensity)))
+        return logPriorDensity;
+    }
+  return logPriorDensity;
+};
+
+/**
+ * Return random number for prior fromits distribution
+ */
+void
+LogPriorDensity::computeNewParams(Vector &ep)
+{
+
+}

mex/sources/estimation/LogPriorDensity.hh

+///////////////////////////////////////////////////////////
+//  LogPriorDensity.hh
+//  Implementation of the Class LogPriorDensity
+//  Created on:      10-Feb-2010 20:56:08
+///////////////////////////////////////////////////////////
+
+#if !defined(LPD_011FD4CF_17CE_4805_882B_046AA07CF443__INCLUDED_)
+#define LPD_011FD4CF_17CE_4805_882B_046AA07CF443__INCLUDED_
+
+//#include <boost/random/variate_generator.hpp>
+#include "Vector.hh"
+#include "EstimatedParametersDescription.hh"
+
+class LogPriorDensity {
+
+public:
+  LogPriorDensity(EstimatedParametersDescription& estParsDesc);
+  virtual ~LogPriorDensity();
+
+  double compute(const Vector &estParams);
+  double
+  getLogPriorDensity()
+  {
+    return logPriorDensity;
+  };
+  void computeNewParams(Vector &newParams);
+
+private:
+  double logPriorDensity;
+  const EstimatedParametersDescription &estParsDesc;
+
+};
+
+#endif // !defined(011FD4CF_17CE_4805_882B_046AA07CF443__INCLUDED_)

mex/sources/estimation/Prior.cc

@@ -25,8 +25,8 @@
 
 #include "Prior.hh"
 
-Prior::Prior(Prior::pShape shape_arg, double mean_arg, double mode_arg, double standard_arg, double lower_bound_arg,      double upper_bound_arg, double fhp_arg, double shp_arg) :
-  shape(shape_arg), mean(mean_arg), mode(mode_arg), standard(standard_arg), lower_bound(lower_bound_arg),       upper_bound(upper_bound_arg),   fhp(fhp_arg), shp(shp_arg)
+Prior::Prior(double mean_arg, double mode_arg, double standard_arg, double lower_bound_arg,      double upper_bound_arg, double fhp_arg, double shp_arg) :
+  mean(mean_arg), mode(mode_arg), standard(standard_arg), lower_bound(lower_bound_arg),       upper_bound(upper_bound_arg),   fhp(fhp_arg), shp(shp_arg)
 {
 
 }

mex/sources/estimation/Prior.hh

 /*
- * Copyright (C) 2009-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/>.
- */
+* Copyright (C) 2009-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/>.
+*/
 
 ///////////////////////////////////////////////////////////
 //  Prior.h
@@ -26,39 +26,208 @@
 #if !defined(Prior_8D5F562F_C831_43f3_B390_5C4EF4433756__INCLUDED_)
 #define Prior_8D5F562F_C831_43f3_B390_5C4EF4433756__INCLUDED_
 
+#include <boost/math/distributions/beta.hpp> // for beta_distribution.
+#include <boost/math/distributions/gamma.hpp> // for gamma_distribution.
+#include <boost/math/distributions/normal.hpp> // for normal_distribution.
+#include <boost/math/distributions/uniform.hpp> // for uniform_distribution.
+
+
 struct Prior
-{
+  {
+
+  public:
+    //! probablity density functions
+    enum pShape
+      {
+      Beta = 1,
+      Gamma = 2,
+      Gaussian = 3, // i.e. Normal density
+      Inv_gamma_1 = 4, // Inverse gamma (type 1) density
+      Uniform = 5,
+      Inv_gamma_2 = 6 //Inverse gamma (type 2) density
+      };
+
+    Prior(double mean, double mode, double standard, double lower_bound, double upper_bound, double fhp, double shp);
+    virtual ~Prior();
+
+    const double mean;
+    const double mode;
+    const double standard;
+    const double lower_bound;
+    const double upper_bound;
+    /**
+    * first  shape parameter
+    */
+    const double fhp;
+    /**
+    * second shape parameter
+    */
+    const double shp;
+    virtual pShape getShape()=0; // e.g. = Beta for  beta shape?
+
+    virtual double pdf(double x) // probability density function value for x
+      {
+      std::cout << "Parent pdf undefined at parent level" << std::endl ;
+      return 0.0;
+      };
+
+    virtual double drand(double x) // rand for density 
+      {
+      std::cout << "Parent rand undefined at parent level" << std::endl ;
+      return 0.0;
+      };
+
+  };
+
+struct BetaPrior: public Prior
+  {
+  public:
+    boost::math::beta_distribution<double> distribution;
+
+    BetaPrior( double mean, double mode, double standard, double lower_bound, double upper_bound, double fhp, double shp)
+      : Prior( mean, mode, standard, lower_bound, upper_bound, fhp, shp),
+      distribution(fhp, shp)
+      {};
+    virtual ~BetaPrior(){};
+    virtual pShape getShape(){return Prior::Beta;}; // e.g. = Beta for  beta shape?
 
-public:
-  //! probablity density functions
-  enum pShape
+    virtual double pdf(double x)
+      {
+      double scalled=x;
+      if (lower_bound ||1.0-upper_bound)
+        scalled=(x- lower_bound)/(upper_bound- lower_bound);    
+      return boost::math::pdf(distribution,scalled);
+      };
+
+    virtual double drand(double x) // rand for density 
+      {
+      return 0.0;
+      };
+  };
+
+struct GammaPrior: public Prior
   {
-    Beta = 1,
-    Gamma = 2,
-    Gaussian = 3, // i.e. Normal density
-    Inv_gamma_1 = 4, // Inverse gamma (type 1) density
-    Uniform = 5,
-    Inv_gamma_2 = 6 //Inverse gamma (type 2) density
+  public:
+    boost::math::gamma_distribution<double> distribution;
+
+    GammaPrior( double mean, double mode, double standard, 
+      double lower_bound, double upper_bound, double fhp, double shp)
+      : Prior( mean, mode, standard, lower_bound, upper_bound, fhp, shp),
+      distribution(fhp, shp)
+      {};
+    virtual ~GammaPrior(){};
+    virtual pShape getShape(){return Prior::Gamma;}; // e.g. = Beta for  beta shape?
+    virtual   double pdf(double x)
+      {
+      return boost::math::pdf(distribution,x- lower_bound);
+      };
+    virtual double drand(double x) // rand for density 
+      {
+      return 0.0;
+      };
   };
 
-  Prior(Prior::pShape shape, double mean, double mode, double standard, double lower_bound, double upper_bound,  double fhp,     double shp);
-  virtual ~Prior();
-
-  const pShape shape;
-  const double mean;
-  const double mode;
-  const double standard;
-  const double lower_bound;
-  const double upper_bound;
-  /**
-   * first  shape parameter
-   */
-  const double fhp;
-  /**
-   * second shape parameter
-   */
-  const double shp;
-
-};
+//  X ~ IG1(s,nu) if X = sqrt(Y) where Y ~ IG2(s,nu) and Y = inv(Z) with Z ~ G(nu/2,2/s) (Gamma distribution) 
+// lpdfig1(x,s,n)= lpdfgam(1/(x*x),n/2,2/s)-2*log(x*x)
+struct InvGamma1_Prior: public Prior
+  {
+  public:
+    boost::math::gamma_distribution<double> distribution;
+    InvGamma1_Prior( double mean, double mode, double standard, 
+      double lower_bound, double upper_bound, double fhp, double shp)
+      : Prior( mean, mode, standard, lower_bound, upper_bound, fhp, shp),
+      distribution(shp/2, 2/fhp)
+      {};
+    virtual ~InvGamma1_Prior(){};
+    virtual pShape getShape(){return Prior::Inv_gamma_1;}; // e.g. = Beta for  beta shape?
+    virtual  double pdf(double x)
+      {
+      double scalled=((x- lower_bound)*(x-lower_bound));
+      if (x>lower_bound)
+        return boost::math::pdf(distribution,1/scalled) / (scalled*scalled);
+      else
+        return 0;
+      };
+    virtual double drand(double x) // rand for density 
+      {
+      return 0.0;
+      };
+  };
+
+// If x~InvGamma(a,b) , then  1/x ~Gamma(a,1/b) distribution
+struct InvGamma2_Prior: public Prior
+  {
+  public:
+    boost::math::gamma_distribution<double> distribution;
+
+    InvGamma2_Prior( double mean, double mode, double standard, 
+      double lower_bound, double upper_bound, double fhp, double shp)
+      : Prior( mean, mode, standard, lower_bound, upper_bound, fhp, shp),
+      distribution(shp/2, 2/fhp)
+      {};
+    virtual ~InvGamma2_Prior(){};
+    virtual pShape getShape(){return Prior::Inv_gamma_2;}; // e.g. = Beta for  beta shape?
+
+    virtual   double pdf(double x)
+      {
+      double scalled= x - lower_bound;
+      if (scalled>0)
+        return boost::math::pdf(distribution,1/scalled) / (scalled*scalled);
+      else
+        return 0;
+      };
+    virtual double drand(double x) // rand for density 
+      {
+      return 0.0;
+      };
+  };
+
+struct GaussianPrior: public Prior
+  {
+  public:
+    boost::math::normal_distribution<double> distribution;
+
+    GaussianPrior( double mean, double mode, double standard, double lower_bound, double upper_bound, double fhp, double shp)
+      : Prior( mean, mode, standard, lower_bound, upper_bound, fhp, shp),
+      distribution(fhp, shp) //distribution(mean, standard)
+      {};
+    virtual ~GaussianPrior(){};
+    virtual pShape getShape(){return Prior::Gaussian;}; // e.g. = Beta for  beta shape?
+    virtual  double pdf(double x)
+      {
+      if (x>lower_bound && x<upper_bound )
+        return boost::math::pdf(distribution,x);
+      else
+        return 0;
+      };
+    virtual double drand(double x) // rand for density 
+      {
+      return 0.0;
+      };
+  };
+
+struct UniformPrior: public Prior
+  {
+  public:
+    boost::math::uniform_distribution<double> distribution;
+
+    UniformPrior( double mean, double mode, double standard, double lower_bound, double upper_bound, double fhp, double shp)
+      : Prior( mean, mode, standard, lower_bound, upper_bound, fhp, shp),
+      distribution(fhp, shp) //distribution(lower_bound, upper_bound)
+      {};
+    virtual ~UniformPrior(){};
+    virtual pShape getShape(){return Prior::Uniform;}; // e.g. = Beta for  beta shape?
+    virtual   double pdf(double x)
+      {
+      if (x>lower_bound && x<upper_bound )
+        return boost::math::pdf(distribution,x);
+      else
+        return 0;
+      };
+    virtual double drand(double x) // rand for density 
+      {
+      return 0.0;
+      };
+  };
 
 #endif // !defined(8D5F562F_C831_43f3_B390_5C4EF4433756__INCLUDED_)