C++ Estimation DLL: Update of logMHMCMCposterior.cc with a draft octave MAT...

C++ Estimation DLL: Update of logMHMCMCposterior.cc with a draft octave MAT draws file save and adding a test random_walk_metropolis_hastings_core.m: Octave version crashes at start of DLL and Matlab version finishes with low acceptance due to frequent B&K and reports error within debugger too - needs more debugging!

mex/sources/estimation/KalmanFilter.cc

@@ -80,15 +80,24 @@ KalmanFilter::compute(const MatrixConstView &dataView, VectorView &steadyState,
                       VectorView &vll, MatrixView &detrendedDataView,
                       size_t start, size_t period, double &penalty, int &info)
 {
-  if(period==0) // initialise all KF matrices
-    initKalmanFilter.initialize(steadyState, deepParams, R, Q, RQRt, T, Pstar, Pinf,
-                              penalty, dataView, detrendedDataView, info);
-  else  // initialise parameter dependent KF matrices only but not Ps
-    initKalmanFilter.initialize(steadyState, deepParams, R, Q, RQRt, T, 
-                              penalty, dataView, detrendedDataView, info);
-
-  double lik= filter(detrendedDataView, H, vll, start, info);
-
+  double lik=INFINITY;
+  try
+    {
+      if(period==0) // initialise all KF matrices
+        initKalmanFilter.initialize(steadyState, deepParams, R, Q, RQRt, T, Pstar, Pinf,
+                                  penalty, dataView, detrendedDataView, info);
+      else  // initialise parameter dependent KF matrices only but not Ps
+        initKalmanFilter.initialize(steadyState, deepParams, R, Q, RQRt, T, 
+                                  penalty, dataView, detrendedDataView, info);
+
+      lik= filter(detrendedDataView, H, vll, start, info);
+    }
+  catch (const DecisionRules::BlanchardKahnException &bke)
+    {
+      info =22;
+      return penalty;
+    }
+ 
   if (info != 0)
     return penalty;
   else 
@@ -97,7 +106,7 @@ KalmanFilter::compute(const MatrixConstView &dataView, VectorView &steadyState,
 };
 
 /**
- * 30:*
+ * Multi-variate standard Kalman Filter
  */
 double
 KalmanFilter::filter(const MatrixView &detrendedDataView,  const Matrix &H, VectorView &vll, size_t start, int &info)
@@ -110,6 +119,7 @@ KalmanFilter::filter(const MatrixView &detrendedDataView,  const Matrix &H, Vect
       if (nonstationary)
         {
           // K=PZ'
+          //blas::gemm("N", "T", 1.0, Pstar, Z, 0.0, K);
           blas::symm("L", "U", 1.0, Pstar, Zt, 0.0, K);
 
           //F=ZPZ' +H = ZK+H
@@ -141,6 +151,7 @@ KalmanFilter::filter(const MatrixView &detrendedDataView,  const Matrix &H, Vect
                   Pstar(i,j)*=0.5;
 
               // K=PZ'
+              //blas::gemm("N", "T", 1.0, Pstar, Z, 0.0, K);
               blas::symm("L", "U", 1.0, Pstar, Zt, 0.0, K);
 
               //F=ZPZ' +H = ZK+H
@@ -167,7 +178,8 @@ KalmanFilter::filter(const MatrixView &detrendedDataView,  const Matrix &H, Vect
           Fdet = 1;
           for (size_t d = 1; d <= p; ++d)
             Fdet *= FUTP(d + (d-1)*d/2 -1);
-          Fdet *=Fdet;
+          Fdet *=Fdet;//*pow(-1.0,p);
+
           logFdet=log(fabs(Fdet));
 
           Ptmp = Pstar;
@@ -181,10 +193,6 @@ KalmanFilter::filter(const MatrixView &detrendedDataView,  const Matrix &H, Vect
           // 3) Pt+1= Ptmp*T' +RQR'
           Pstar = RQRt;
           blas::gemm("N", "T", 1.0, Ptmp, T, 1.0, Pstar);
-          //enforce Pstar symmetry with P=(P+P')/2=0.5P+0.5P'
-          //blas::gemm("N", "T", 0.5, Ptmp, T, 0.5, Pstar);
-          //mat::transpose(Ptmp, Pstar);
-          //mat::add(Pstar,Ptmp);
 
           if (t>0)
             nonstationary = mat::isDiff(KFinv, oldKFinv, riccati_tol);

mex/sources/estimation/LogLikelihoodSubSample.cc

@@ -36,7 +36,7 @@ LogLikelihoodSubSample::LogLikelihoodSubSample(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 double qz_criterium,
                                                const std::vector<size_t> &varobs, double riccati_tol, double lyapunov_tol, int &INinfo) :
-  estiParDesc(INestiParDesc),
+  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)
 {
@@ -46,6 +46,8 @@ double
 LogLikelihoodSubSample::compute(VectorView &steadyState, const MatrixConstView &dataView, const Vector &estParams, Vector &deepParams,
                                 Matrix &Q, Matrix &H, VectorView &vll, MatrixView &detrendedDataView, int &info, size_t start, size_t period)
 {
+  penalty=startPenalty;
+  logLikelihood=startPenalty;
 
   updateParams(estParams, deepParams, Q, H, period);
   if (info == 0)
@@ -84,17 +86,11 @@ LogLikelihoodSubSample::updateParams(const Vector &estParams, Vector &deepParams
               break;
 
             case EstimatedParameter::measureErr_SD:
-#ifdef DEBUG
-              mexPrintf("Setting of H var_endo\n");
-#endif
               k = estiParDesc.estParams[i].ID1;
               H(k, k) = estParams(i)*estParams(i);
               break;
 
             case EstimatedParameter::shock_Corr:
-#ifdef DEBUG
-              mexPrintf("Setting of Q corrx\n");
-#endif
               k1 = estiParDesc.estParams[i].ID1;
               k2 = estiParDesc.estParams[i].ID2;
               Q(k1, k2) = estParams(i)*sqrt(Q(k1, k1)*Q(k2, k2));
@@ -128,9 +124,6 @@ LogLikelihoodSubSample::updateParams(const Vector &estParams, Vector &deepParams
               break;
 
             case EstimatedParameter::measureErr_Corr:
-#ifdef DEBUG
-              mexPrintf("Setting of H corrn\n");
-#endif
               k1 = estiParDesc.estParams[i].ID1;
               k2 = estiParDesc.estParams[i].ID2;
               //      H(k1,k2) = xparam1(i)*sqrt(H(k1,k1)*H(k2,k2));

mex/sources/estimation/LogLikelihoodSubSample.hh

@@ -43,13 +43,14 @@ public:
   virtual ~LogLikelihoodSubSample();
 
 private:
-  double penalty;
+  double startPenalty, penalty;
   double logLikelihood;
   EstimatedParametersDescription &estiParDesc;
   KalmanFilter kalmanFilter;
   VDVEigDecomposition eigQ;
   VDVEigDecomposition eigH;
   int &info;
+
   // methods
   void updateParams(const Vector &estParams, Vector &deepParams,
                     Matrix &Q, Matrix &H, size_t period);

mex/sources/estimation/RandomWalkMetropolisHastings.cc

@@ -19,17 +19,27 @@
 
 #include "RandomWalkMetropolisHastings.hh"
 
+#include <iostream>
+#include <fstream>
+
 double
 RandomWalkMetropolisHastings::compute(VectorView &mhLogPostDens, MatrixView &mhParams, Matrix &steadyState,
                                       Vector &estParams, Vector &deepParams, const MatrixConstView &data, Matrix &Q, Matrix &H,
                                       const size_t presampleStart, int &info, const size_t startDraw, size_t nMHruns, const Matrix &Dscale,
                                       LogPosteriorDensity &lpd, Prior &drawDistribution, EstimatedParametersDescription &epd)
 {
+  //streambuf *likbuf, *drawbuf *backup;
+  std::ofstream urandfilestr, drawfilestr;
+  urandfilestr.open ("urand.csv");
+  drawfilestr.open ("paramdraws.csv");
+
   bool overbound;
-  double newLogpost, logpost;
+  double newLogpost, logpost, urand;
   size_t count, accepted = 0;
   parDraw = estParams;
+
   logpost =  - lpd.compute(steadyState, estParams, deepParams, data, Q, H, presampleStart, info);
+
   for (size_t run = startDraw - 1; run < nMHruns; ++run)
     {
       overbound=false;
@@ -49,25 +59,36 @@ RandomWalkMetropolisHastings::compute(VectorView &mhLogPostDens, MatrixView &mhP
             {
               newLogpost = - lpd.compute(steadyState, newParDraw, deepParams, data, Q, H, presampleStart, info);
             }
-          catch(...)
+          catch(const std::exception &e)
+            {
+              throw; // for now handle the system and other errors higher-up
+            }
+          catch (...)
             {
               newLogpost = -INFINITY;
             }
         }
-      if ((newLogpost > -INFINITY) && log(uniform.drand()) < newLogpost-logpost)
+      urand=uniform.drand();
+      if ((newLogpost > -INFINITY) && log(urand) < newLogpost-logpost)
         {
-          mat::get_row(mhParams, run) = newParDraw;
           parDraw = newParDraw;
-          mhLogPostDens(run) = newLogpost;
           logpost = newLogpost;
           accepted++;
         }
-      else
-        {
-          mat::get_row(mhParams, run) = parDraw;
-          mhLogPostDens(run) = logpost;
-        }
+      mat::get_row(mhParams, run) = parDraw;
+      mhLogPostDens(run) = logpost;
+
+      //urandfilestr.write(urand);
+      urandfilestr << urand << "\n"; //","
+      for (size_t c=0;c<newParDraw.getSize()-1;++c)
+        drawfilestr << newParDraw(c) << ",";
+//        drawfilestr.write(newParDraw(i));
+      drawfilestr <<  newParDraw(newParDraw.getSize()-1) << "\n";
     }
+
+  urandfilestr.close();
+  drawfilestr.close();
+
   return (double) accepted/(nMHruns-startDraw+1);
 }