diff --git a/mex/sources/estimation/InitializeKalmanFilter.cc b/mex/sources/estimation/InitializeKalmanFilter.cc
index 67eb8ef8ce69333ae4a80cf586584a0f1c180004..76e5d7330c79c6d43d9c16bfc1816caff7eb116e 100644
--- a/mex/sources/estimation/InitializeKalmanFilter.cc
+++ b/mex/sources/estimation/InitializeKalmanFilter.cc
@@ -50,10 +50,10 @@ InitializeKalmanFilter::InitializeKalmanFilter(const std::string &modName, size_
if (inv_order_var_arg[i]-n_static < n_pred && inv_order_var_arg[i]-n_static >= 0)
inv_ric[j++] = ricIN[inv_order_var_arg[i]-n_static];
}
-
+// initialise parameter dependent KF matrices only but not Ps
void
InitializeKalmanFilter::initialize(Vector &steadyState, const Vector &deepParams, Matrix &R, const Matrix &Z,
- const Matrix &Q, Matrix &RQRt, Matrix &T, Matrix &Pstar, Matrix &Pinf,
+ const Matrix &Q, Matrix &RQRt, Matrix &T,
double &penalty, const MatrixView &dataView, Matrix &Y, int &info)
{
modelSolution.compute(steadyState, deepParams, ghx_raw, ghu_raw);
@@ -63,10 +63,20 @@ InitializeKalmanFilter::initialize(Vector &steadyState, const Vector &deepParams
mat::reorderRowsByVectors(ghu, mat::nullVec, ghu_raw, order_var);
setT(T, info);
- setR(R, info);
- setPstar(Pstar, Pinf, T, R, Q, RQRt, info);
+ setRQR(R, Q, RQRt, info);
+}
+
+// initialise all KF matrices
+void
+InitializeKalmanFilter::initialize(Vector &steadyState, const Vector &deepParams, Matrix &R, const Matrix &Z,
+ const Matrix &Q, Matrix &RQRt, Matrix &T, Matrix &Pstar, Matrix &Pinf,
+ double &penalty, const MatrixView &dataView, Matrix &Y, int &info)
+{
+ initialize(steadyState, deepParams, R, Z, Q, RQRt, T, penalty, dataView, Y, info);
+ setPstar(Pstar, Pinf, T, RQRt, info);
}
+
void
InitializeKalmanFilter::setT(Matrix &T, int &info)
{
@@ -79,22 +89,22 @@ InitializeKalmanFilter::setT(Matrix &T, int &info)
}
void
-InitializeKalmanFilter::setR(Matrix &R, int &info)
+InitializeKalmanFilter::setRQR(Matrix &R, const Matrix &Q, Matrix &RQRt, int &info)
{
R.setAll(0.0);
//B(i_n_iv,:) = dr.ghu(iv,:); and R=B;
mat::assignByVectors(R, mat::nullVec, mat::nullVec, ghu, riv, mat::nullVec);
-}
-
-void
-InitializeKalmanFilter::setPstar(Matrix &Pstar, Matrix &Pinf, Matrix &T, Matrix &R, const Matrix &Q, Matrix &RQRt, int &info)
-{
// Matrix RQRt=R*Q*R'
RQ.setAll(0.0);
blas::gemm("N", "N", 1.0, R, Q, 1.0, RQ); // R*Q
RQRt.setAll(0.0);
blas::gemm("N", "T", 1.0, RQ, R, 1.0, RQRt); // R*Q*R'
+}
+
+void
+InitializeKalmanFilter::setPstar(Matrix &Pstar, Matrix &Pinf, const Matrix &T, const Matrix &RQRt, int &info)
+{
try
{
diff --git a/mex/sources/estimation/InitializeKalmanFilter.hh b/mex/sources/estimation/InitializeKalmanFilter.hh
index fac1ed90eff0f04985978794e2c0fe5e54d1e6e6..ab227b6c12f4d74e3d12d2b5cd7a4a45c80f4f80 100644
--- a/mex/sources/estimation/InitializeKalmanFilter.hh
+++ b/mex/sources/estimation/InitializeKalmanFilter.hh
@@ -45,8 +45,12 @@ public:
const Matrix &llincidence, double qz_criterium, const std::vector<size_t> &order_var_arg, const std::vector<size_t> &inv_order_var_arg,
const std::vector<size_t> &riv, const std::vector<size_t> &ric, double lyapunov_tol, int &info);
virtual ~InitializeKalmanFilter();
+ // initialise all KF matrices
void initialize(Vector &steadyState, const Vector &deepParams, Matrix &R, const Matrix &Z, const Matrix &Q, Matrix &RQRt,
Matrix &T, Matrix &Pstar, Matrix &Pinf, double &penalty, const MatrixView &dataView, Matrix &Y, int &info);
+ // initialise parameter dependent KF matrices only but not Ps
+ void initialize(Vector &steadyState, const Vector &deepParams, Matrix &R, const Matrix &Z, const Matrix &Q, Matrix &RQRt,
+ Matrix &T, double &penalty, const MatrixView &dataView, Matrix &Y, int &info);
private:
const std::vector<size_t> riv; // restrict_var_list
@@ -61,8 +65,8 @@ private:
Matrix ghu, ghu_raw;
Matrix Rt, RQ;
void setT(Matrix &T, int &info);
- void setR(Matrix &R, int &info);
- void setPstar(Matrix &Pstar, Matrix &Pinf, Matrix &T, Matrix &R, const Matrix &Q, Matrix &RQRt, int &info);
+ void setRQR(Matrix &R, const Matrix &Q, Matrix &RQRt, int &info);
+ void setPstar(Matrix &Pstar, Matrix &Pinf, const Matrix &T, const Matrix &RQRt, int &info);
};
diff --git a/mex/sources/estimation/KalmanFilter.cc b/mex/sources/estimation/KalmanFilter.cc
index e208ccd7d589dddb9583e4676ebe54a502ad7bfc..9b4510f809d8c5e67e2e4a14c13b8524e16050b0 100644
--- a/mex/sources/estimation/KalmanFilter.cc
+++ b/mex/sources/estimation/KalmanFilter.cc
@@ -56,12 +56,16 @@ KalmanFilter::KalmanFilter(const std::string &modName, size_t n_endo, size_t n_e
double
KalmanFilter::compute(const MatrixView &dataView, Vector &steadyState,
const Matrix &Q, const Matrix &H, const Vector &deepParams,
- Vector &vll, size_t start, double &penalty, int &info)
+ VectorView &vll, size_t start, size_t period, double &penalty, int &info)
{
double ll;
Matrix Y(dataView.getRows(), dataView.getCols()); // data
- initKalmanFilter.initialize(steadyState, deepParams, R, Z, Q, RQRt, T, Pstar, Pinf,
+ if(period==0) // initialise all KF matrices
+ initKalmanFilter.initialize(steadyState, deepParams, R, Z, Q, RQRt, T, Pstar, Pinf,
+ penalty, dataView, Y, info);
+ else // initialise parameter dependent KF matrices only but not Ps
+ initKalmanFilter.initialize(steadyState, deepParams, R, Z, Q, RQRt, T,
penalty, dataView, Y, info);
return ll = filter(Y, H, vll, start, info);
@@ -72,7 +76,7 @@ KalmanFilter::compute(const MatrixView &dataView, Vector &steadyState,
* 30:*
*/
double
-KalmanFilter::filter(const Matrix &dataView, const Matrix &H, Vector &vll, size_t start, int &info)
+KalmanFilter::filter(const Matrix &dataView, const Matrix &H, VectorView &vll, size_t start, int &info)
{
double loglik=0.0, ll, logFdet, Fdet;
int p = Finv.getRows();
diff --git a/mex/sources/estimation/KalmanFilter.hh b/mex/sources/estimation/KalmanFilter.hh
index fab350f35e41d42860566a51c6e6e39faeadaa96..c65e665a0a61f5cc976862ae6e073ba83915b322 100644
--- a/mex/sources/estimation/KalmanFilter.hh
+++ b/mex/sources/estimation/KalmanFilter.hh
@@ -56,7 +56,7 @@ public:
double compute(const MatrixView &dataView, Vector &steadyState,
const Matrix &Q, const Matrix &H, const Vector &deepParams,
- Vector &vll, size_t start, double &penalty, int &info);
+ VectorView &vll, size_t start, size_t period, double &penalty, int &info);
protected:
// static double calcStepLogLik(const PLUFact &Finv, const Vector &v);
@@ -77,7 +77,7 @@ private:
Matrix vtFinv, vtFinvVt; // intermeiate observation error *Finv vector
double riccati_tol;
InitializeKalmanFilter initKalmanFilter; //Initialise KF matrices
- double filter(const Matrix &dataView, const Matrix &H, Vector &vll, size_t start, int &info);
+ double filter(const Matrix &dataView, const Matrix &H, VectorView &vll, size_t start, int &info);
};
diff --git a/mex/sources/estimation/LogLikelihoodSubSample.cc b/mex/sources/estimation/LogLikelihoodSubSample.cc
new file mode 100644
index 0000000000000000000000000000000000000000..ab66589db23489885d5bc837686aeccf4297a43a
--- /dev/null
+++ b/mex/sources/estimation/LogLikelihoodSubSample.cc
@@ -0,0 +1,50 @@
+/*
+ * 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/>.
+ */
+
+///////////////////////////////////////////////////////////
+// LogLikelihoodSubSample.cpp
+// Implementation of the Class LogLikelihoodSubSample
+// Created on: 14-Jan-2010 22:39:14
+///////////////////////////////////////////////////////////
+
+#include "LogLikelihoodSubSample.hh"
+
+LogLikelihoodSubSample::~LogLikelihoodSubSample()
+{
+};
+
+LogLikelihoodSubSample::LogLikelihoodSubSample(const std::string &modName, 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 Matrix &ll_incidence, const double qz_criterium, const std::vector<size_t> &order_var,
+ const std::vector<size_t> &inv_order_var, const std::vector<size_t> &varobs, const std::vector<size_t> &riv,
+ const std::vector<size_t> &ric, double riccati_tol, double lyapunov_tol, int &info) :
+ kalmanFilter(modName, n_endo, n_exo, zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg, zeta_static_arg, ll_incidence, qz_criterium,
+ order_var, inv_order_var, varobs, riv, ric, riccati_tol, lyapunov_tol, info)
+{
+};
+
+double
+LogLikelihoodSubSample::compute(Vector &steadyState, const MatrixView &dataView, const Vector &deepParams, //const estPeriod &estPeriod,
+ VectorView &vll, int &info, size_t start, size_t period, const Matrix &Q, const Matrix &H)
+{
+ logLikelihood = kalmanFilter.compute(dataView, steadyState, Q, H, deepParams, vll, start, period, penalty, info);
+ return logLikelihood;
+
+};
+
diff --git a/mex/sources/estimation/LogLikelihoodSubSample.hh b/mex/sources/estimation/LogLikelihoodSubSample.hh
new file mode 100644
index 0000000000000000000000000000000000000000..06d115eab8f40ee039877209709689bbc861ab4f
--- /dev/null
+++ b/mex/sources/estimation/LogLikelihoodSubSample.hh
@@ -0,0 +1,53 @@
+/*
+ * 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/>.
+ */
+
+///////////////////////////////////////////////////////////
+// LogLikelihoodSubSample.h
+// Implementation of the Class LogLikelihoodSubSample
+// Created on: 14-Jan-2010 22:39:14
+///////////////////////////////////////////////////////////
+
+#if !defined(DF8B7AF5_8169_4587_9037_2CD2C82E2DDF__INCLUDED_)
+#define DF8B7AF5_8169_4587_9037_2CD2C82E2DDF__INCLUDED_
+
+#include "KalmanFilter.hh"
+
+class LogLikelihoodSubSample {
+
+public:
+ LogLikelihoodSubSample(const std::string &modName, 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 Matrix &ll_incidence, const double qz_criterium, const std::vector<size_t> &order_var,
+ const std::vector<size_t> &inv_order_var, const std::vector<size_t> &varobs, const std::vector<size_t> &riv,
+ const std::vector<size_t> &ric, double riccati_tol_in, double lyapunov_tol, int &info);
+
+ double compute(Vector &steadyState, const MatrixView &dataView, const Vector &deepParams, //const estPeriod &estPeriod,
+ VectorView &vll, int &info, size_t start, size_t period, const Matrix &Q, const Matrix &H);
+ virtual ~LogLikelihoodSubSample();
+
+public:
+ KalmanFilter kalmanFilter;
+
+private:
+ double penalty;
+ double logLikelihood;
+
+};
+
+#endif // !defined(DF8B7AF5_8169_4587_9037_2CD2C82E2DDF__INCLUDED_)
diff --git a/mex/sources/estimation/tests/testInitKalman.cc b/mex/sources/estimation/tests/testInitKalman.cc
index ef42863cf0b1083f47206f684fffce83ae85f46e..94f6975dde3e3639a3507159cf650b9a62f05e8a 100644
--- a/mex/sources/estimation/tests/testInitKalman.cc
+++ b/mex/sources/estimation/tests/testInitKalman.cc
@@ -44,15 +44,26 @@ main(int argc, char **argv)
steadyState(n_endo), deepParams(npar);
double dYSparams [] = {
- 1.0002, 0.9933, 1.0070, 1.0000,
- 2.7186, 1.0073, 18.9822, 0.8608,
- 0.3167, 0.8610, 1.0085, 0.9917,
- 1.3559, 1.0085, 0.9929
+ 1.000199998312523,
+ 0.993250551764778,
+ 1.006996670195112,
+ 1,
+ 2.718562165733039,
+ 1.007250753636589,
+ 18.982191739915155,
+ 0.860847884886309,
+ 0.316729149714572,
+ 0.861047883198832,
+ 1.00853622757204,
+ 0.991734328394345,
+ 1.355876776121869,
+ 1.00853622757204,
+ 0.992853374047708
};
double vcov[] = {
- 0.0013, 0,
- 0, 0.0001
+ 0.001256631601, 0.0,
+ 0.0, 0.000078535044
};
Matrix
diff --git a/mex/sources/estimation/tests/testKalman.cc b/mex/sources/estimation/tests/testKalman.cc
index b95504e8511dd3f085747bcf997f72295cf473f4..599085d7b2018cad17a25b7c1c8bda94a4bdefa0 100644
--- a/mex/sources/estimation/tests/testKalman.cc
+++ b/mex/sources/estimation/tests/testKalman.cc
@@ -44,15 +44,26 @@ main(int argc, char **argv)
steadyState(n_endo), deepParams(npar);
double dYSparams [] = {
- 1.0002, 0.9933, 1.0070, 1.0000,
- 2.7186, 1.0073, 18.9822, 0.8608,
- 0.3167, 0.8610, 1.0085, 0.9917,
- 1.3559, 1.0085, 0.9929
+ 1.000199998312523,
+ 0.993250551764778,
+ 1.006996670195112,
+ 1,
+ 2.718562165733039,
+ 1.007250753636589,
+ 18.982191739915155,
+ 0.860847884886309,
+ 0.316729149714572,
+ 0.861047883198832,
+ 1.00853622757204,
+ 0.991734328394345,
+ 1.355876776121869,
+ 1.00853622757204,
+ 0.992853374047708
};
double vcov[] = {
- 0.0013, 0,
- 0, 0.0001
+ 0.001256631601, 0.0,
+ 0.0, 0.000078535044
};
Matrix
@@ -162,12 +173,11 @@ main(int argc, char **argv)
double lyapunov_tol = 1e-16;
double riccati_tol = 1e-16;
int info = 0;
- //const MatrixView dataView(&lyapunov_tol, 1, 1, 1); // dummy
- //Matrix yView(dataView.getRows(),dataView.getCols());
Matrix yView(nobs,192); // dummy
yView.setAll(0.2);
const MatrixView dataView(yView, 0, 0, nobs, yView.getCols() ); // dummy
Vector vll(yView.getCols());
+ VectorView vwll(vll,0,vll.getSize());
double penalty = 1e8;
@@ -177,18 +187,10 @@ main(int argc, char **argv)
std::cout << "Initilise KF with Q: " << std::endl << Q << std::endl;
// std::cout << "and Z" << std::endl << Z << std::endl;
- size_t start=0;
+ size_t start=0, period=0;
double ll=kalman.compute(dataView, steadyState, Q, H, deepParams,
- vll, start , penalty, info);
-
- /*
- std::cout << "Matrix T: " << std::endl << T << std::endl;
- std::cout << "Matrix R: " << std::endl << R << std::endl;
- std::cout << "Matrix RQRt: " << std::endl << RQRt << std::endl;
- std::cout << "Matrix Pstar: " << std::endl << Pstar << std::endl;
- std::cout << "Matrix Pinf: " << std::endl << Pinf << std::endl;
- */
- std::cout << "ll: " << std::endl << ll << std::endl;
+ vwll, start, period, penalty, info);
+ std::cout << "ll: " << std::endl << ll << std::endl;
}