Skip to content
GitLab
Commit 43b2aae6 authored by George Perendia's avatar George Perendia
Browse files

C++ Estimation Minor changes to KalmanFilter initialisation and testing... 

C++ Estimation Minor changes to KalmanFilter initialisation and testing routines to improve handling structural changes parameters and test accuracy
parent e9846ae6
Hide whitespace changes
Inline Side-by-side
mex/sources/estimation/InitializeKalmanFilter.cc
View file @ 43b2aae6
......@@ -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
{
......
mex/sources/estimation/InitializeKalmanFilter.hh
View file @ 43b2aae6
......@@ -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);
};
......
mex/sources/estimation/KalmanFilter.cc
View file @ 43b2aae6
......@@ -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();
......
mex/sources/estimation/KalmanFilter.hh
View file @ 43b2aae6
......@@ -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);
};
......
mex/sources/estimation/LogLikelihoodSubSample.cc 0 → 100644
View file @ 43b2aae6
/*
* 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;
};
mex/sources/estimation/LogLikelihoodSubSample.hh 0 → 100644
View file @ 43b2aae6
/*
* 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_)
mex/sources/estimation/tests/testInitKalman.cc
View file @ 43b2aae6
......@@ -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
......
mex/sources/estimation/tests/testKalman.cc
View file @ 43b2aae6
......@@ -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;
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment