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

C++ Estimation: Adding test and fixing a bug for KalmanFilter.cc and making... 

C++ Estimation: Adding test and fixing a bug for KalmanFilter.cc and making some cosmetic changes to InitializeKalmanFilter.cc and ModelSolution.cc
parent 0f06026c
Hide whitespace changes
Inline Side-by-side
mex/sources/estimation/InitializeKalmanFilter.cc
View file @ a489c881
......@@ -92,11 +92,9 @@ InitializeKalmanFilter::setPstar(Matrix &Pstar, Matrix &Pinf, Matrix &T, Matrix
// Matrix RQRt=R*Q*R'
RQ.setAll(0.0);
blas::gemm("N", "N", 1.0, R, Q, 0.0, RQ); // R*Q
blas::gemm("N", "N", 1.0, R, Q, 1.0, RQ); // R*Q
RQRt.setAll(0.0);
//mat::transpose(Rt, R);
blas::gemm("N", "T", 1.0, RQ, R, 0.0, RQRt); // R*Q*R'
//mat::transpose(RQR);
blas::gemm("N", "T", 1.0, RQ, R, 1.0, RQRt); // R*Q*R'
try
{
......
mex/sources/estimation/KalmanFilter.cc
View file @ a489c881
......@@ -74,7 +74,7 @@ KalmanFilter::compute(const MatrixView &dataView, Vector &steadyState,
double
KalmanFilter::filter(const Matrix &dataView, const Matrix &H, Vector &vll, size_t start, int &info)
{
double loglik, ll, logFdet;
double loglik=0.0, ll, logFdet, Fdet;
int p = Finv.getRows();
bool nonstationary = true;
......@@ -98,9 +98,12 @@ KalmanFilter::filter(const Matrix &dataView, const Matrix &H, Vector &vll, size
KFinv.setAll(0.0);
blas::gemm("N", "N", 1.0, K, Finv, 1.0, KFinv);
// deteminant of F:
logFdet = 1;
Fdet = 1;
for (int d = 0; d < p; ++d)
logFdet *= (-F(d, d));
Fdet *= (-F(d, d));
logFdet=log(fabs(Fdet));
Pold = Pstar;
Ptmp = Pstar;
// Pt+1= T(Pt - KFinvK')T' +RQR'
......@@ -115,9 +118,6 @@ KalmanFilter::filter(const Matrix &dataView, const Matrix &H, Vector &vll, size
blas::gemm("N", "T", 1.0, Ptmp, T, 1.0, Pstar);
nonstationary = mat::isDiffSym(Pstar, Pold, riccati_tol);
}
else
{
}
// err= Yt - Za
//VectorView yt(dataView.getData()+t*dataView.getRows(),dataView.getRows(),1); // current observation vector
......@@ -135,20 +135,13 @@ KalmanFilter::filter(const Matrix &dataView, const Matrix &H, Vector &vll, size
blas::gemm("T", "N", 1.0, vt, Finv, 0.0, vtFinv);
blas::gemm("N", "N", 1.0, vtFinv, vt, 0.0, vtFinvVt);
ll = -0.5*(p*log(2*M_PI)+logFdet+*(vtFinvVt.getData()));
ll = -0.5*(p*log(2*M_PI)+logFdet+(*(vtFinvVt.getData())));
vll(t) = ll;
if (t > start) loglik += ll;
if (t >= start) loglik += ll;
}
return loglik;
}
/**
* 89:*
double KalmanFilter::calcStepLogLik(const PLUFact &Finv, const Vector &v){
}
*/
mex/sources/estimation/ModelSolution.cc
View file @ a489c881
......@@ -62,13 +62,13 @@ ModelSolution::ComputeModelSolution(Vector& steadyState, const Vector& deepParam
Vector llXsteadyState(n_jcols-n_exo);
for (int ll_row = 0; ll_row < ll_incidence.getRows(); ll_row++)
for (size_t ll_row = 0; ll_row < ll_incidence.getRows(); ll_row++)
{
// populate (non-sparse) vector with ysteady values
for (int i = 0; i < n_endo; i++)
for (size_t i = 0; i < n_endo; i++)
{
if (ll_incidence(ll_row, i))
llXsteadyState(((int) ll_incidence(ll_row, i))-1) = steadyState(i);
llXsteadyState(((size_t) ll_incidence(ll_row, i))-1) = steadyState(i);
}
}
#ifdef DEBUG
......@@ -78,8 +78,8 @@ ModelSolution::ComputeModelSolution(Vector& steadyState, const Vector& deepParam
//get jacobian
dynamicDLLp.eval(llXsteadyState, Mx, &deepParams, 1, residual, &jacobian, NULL, NULL);
std::cout << "jacobian: " << std::endl << jacobian << std::endl;
#ifdef DEBUG
std::cout << "jacobian: " << std::endl << jacobian << std::endl;
mexPrintf(" compute rules \n");
#endif
//compute rules
......
mex/sources/estimation/tests/testKalman.cc 0 → 100644
View file @ a489c881
/*
* Copyright (C) 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/>.
*/
// Test for InitializeKalmanFilter
// Uses fs2000k2e.mod and its ..._dynamic.mexw32
#include "KalmanFilter.hh"
int
main(int argc, char **argv)
{
if (argc < 2)
{
std::cerr << argv[0] << ": please provide as argument the name of the dynamic DLL generated from fs2000k2.mod (typically fs2000k2_dynamic.mex*)" << std::endl;
exit(EXIT_FAILURE);
}
std::string modName = argv[1];
const int npar = 7; //(int)mxGetM(mxFldp);
const size_t n_endo = 15, n_exo = 2;
std::vector<size_t> zeta_fwrd_arg;
std::vector<size_t> zeta_back_arg;
std::vector<size_t> zeta_mixed_arg;
std::vector<size_t> zeta_static_arg;
//std::vector<size_t>
double qz_criterium = 1.000001; //1.0+1.0e-9;
Vector
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
};
double vcov1[] = {
0.1960e-3, 0.0,
0.0, 0.0250e-3
};
double vcov[] = {
0.0013, 0,
0, 0.0001
};
Matrix
ll_incidence(3, n_endo); // leads and lags indices
double inllincidence[] = {
1, 5, 0,
2, 6, 20,
0, 7, 21,
0, 8, 0,
0, 9, 0,
0, 10, 0,
3, 11, 0,
0, 12, 0,
0, 13, 0,
0, 14, 0,
0, 15, 0,
0, 16, 0,
4, 17, 0,
0, 18, 0,
0, 19, 22,
};
MatrixView
llincidence(inllincidence, 3, n_endo, 3); // leads and lags indices
ll_incidence = llincidence;
double dparams[] = {
0.3560,
0.9930,
0.0085,
1.0002,
0.1290,
0.6500,
0.0100
};
VectorView
modParamsVW(dparams, npar, 1);
deepParams = modParamsVW;
VectorView
steadyStateVW(dYSparams, n_endo, 1);
steadyState = steadyStateVW;
std::cout << "Vector deepParams: " << std::endl << deepParams << std::endl;
std::cout << "MatrixVw llincidence: " << std::endl << llincidence << std::endl;
std::cout << "Matrix ll_incidence: " << std::endl << ll_incidence << std::endl;
std::cout << "Vector steadyState: " << std::endl << steadyState << std::endl;
// Set zeta vectors [0:(n-1)] from Matlab indices [1:n] so that:
// order_var = [ stat_var(:); pred_var(:); both_var(:); fwrd_var(:)];
size_t statc[] = { 4, 5, 6, 8, 9, 10, 11, 12, 14};
size_t back[] = {1, 7, 13};
size_t both[] = {2};
size_t fwd[] = { 3, 15};
for (int i = 0; i < 9; ++i)
zeta_static_arg.push_back(statc[i]-1);
for (int i = 0; i < 3; ++i)
zeta_back_arg.push_back(back[i]-1);
for (int i = 0; i < 1; ++i)
zeta_mixed_arg.push_back(both[i]-1);
for (int i = 0; i < 2; ++i)
zeta_fwrd_arg.push_back(fwd[i]-1);
size_t nriv = 6, nric = 4;
size_t sriv[] = {7, 8, 10, 11, 12, 13};
size_t sric[] = {3, 4, 5, 6};
std::vector<size_t> riv;
for (size_t i = 0; i < nriv; ++i)
riv.push_back(sriv[i]-1);
std::vector<size_t> ric;
for (size_t i = 0; i < nric; ++i)
ric.push_back(sric[i]-1);
size_t nobs = 2;
size_t svarobs[] = {2, 1};
std::vector<size_t> varobs;
for (size_t i = 0; i < nobs; ++i)
varobs.push_back(svarobs[i]-1);
Matrix Q(n_exo), H(nobs);
/*
T(riv.size(), riv.size()), R(riv.size(), n_exo),
RQRt(riv.size(), riv.size()), Pstar(riv.size(), riv.size()),
Pinf(riv.size(), riv.size()), Z(varobs.size(), riv.size()),
*/
H.setAll(0.0);
/* for (size_t i = 0; i < varobs.size(); ++i)
Z(i, varobs[i]) = 1.0;
*/
MatrixView
vCovVW(vcov, n_exo, n_exo, n_exo);
Q = vCovVW;
std::cout << "Matrix Q: " << std::endl << Q << std::endl;
size_t sorder_var[] =
{ 4, 5, 6, 8, 9, 10, 11, 12, 14, 1, 7, 13, 2, 3, 15};
std::vector<size_t> order_var;
for (size_t ii = 0; ii < n_endo; ++ii)
order_var.push_back(sorder_var[ii]-1); //= (1:endo_nbr)';
size_t sinv_order_var[] =
{ 10, 13, 14, 1, 2, 3, 11, 4, 5, 6, 7, 8, 12, 9, 15};
std::vector<size_t> inv_order_var;
for (size_t ii = 0; ii < n_endo; ++ii)
inv_order_var.push_back(sinv_order_var[ii]-1); //= (1:endo_nbr)';
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());
double penalty = 1e8;
KalmanFilter kalman(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);
std::cout << "Initilise KF with Q: " << std::endl << Q << std::endl;
// std::cout << "and Z" << std::endl << Z << std::endl;
size_t start=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;
}
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