diff --git a/mex/sources/estimation/InitializeKalmanFilter.cc b/mex/sources/estimation/InitializeKalmanFilter.cc index da16e46f5..4d237e450 100644 --- a/mex/sources/estimation/InitializeKalmanFilter.cc +++ b/mex/sources/estimation/InitializeKalmanFilter.cc @@ -33,46 +33,48 @@ InitializeKalmanFilter::~InitializeKalmanFilter() InitializeKalmanFilter::InitializeKalmanFilter(const std::string &modName, size_t n_endo_arg, size_t n_exo_arg, const std::vector &zeta_fwrd_arg, const std::vector &zeta_back_arg, const std::vector &zeta_mixed_arg, const std::vector &zeta_static_arg, - double qz_criterium, const std::vector &order_var_arg, const std::vector &inv_order_var_arg, - const std::vector &rivIN, const std::vector &ricIN, double lyapunov_tolIN, int &info) : - riv(rivIN), inv_ric(ricIN.size()), order_var(order_var_arg), lyapunov_tol(lyapunov_tolIN), + const std::vector &zeta_varobs_back_mixed_arg, + double qz_criterium_arg, + double lyapunov_tol_arg, int &info) : + lyapunov_tol(lyapunov_tol_arg), + zeta_varobs_back_mixed(zeta_varobs_back_mixed_arg), detrendData(false), modelSolution(modName, n_endo_arg, n_exo_arg, zeta_fwrd_arg, zeta_back_arg, - zeta_mixed_arg, zeta_static_arg, qz_criterium), discLyapFast(riv.size()), - ghx(n_endo_arg, zeta_back_arg.size() + zeta_mixed_arg.size()), - ghx_raw(n_endo_arg, zeta_back_arg.size() + zeta_mixed_arg.size()), - ghu(n_endo_arg, n_exo_arg), ghu_raw(n_endo_arg, n_exo_arg), - Rt(n_exo_arg, riv.size()), RQ(riv.size(), n_exo_arg) + zeta_mixed_arg, zeta_static_arg, qz_criterium_arg), + discLyapFast(zeta_varobs_back_mixed.size()), + g_x(n_endo_arg, zeta_back_arg.size() + zeta_mixed_arg.size()), + g_u(n_endo_arg, n_exo_arg), + Rt(n_exo_arg, zeta_varobs_back_mixed.size()), + RQ(zeta_varobs_back_mixed.size(), n_exo_arg) { - size_t n_pred = ghx.getCols(); - size_t n_static = zeta_static_arg.size(); - size_t j = 0; - for (size_t i = 0; i < n_endo_arg; ++i) - 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]; + std::vector zeta_back_mixed; + set_union(zeta_back_arg.begin(), zeta_back_arg.end(), + zeta_mixed_arg.begin(), zeta_mixed_arg.end(), + back_inserter(zeta_back_mixed)); + for (size_t i = 0; i < zeta_back_mixed.size(); i++) + pi_bm_vbm.push_back(find(zeta_varobs_back_mixed.begin(), zeta_varobs_back_mixed.end(), + zeta_back_mixed[i]) - zeta_varobs_back_mixed.begin()); + } // initialise parameter dependent KF matrices only but not Ps void -InitializeKalmanFilter::initialize(Vector &steadyState, const Vector &deepParams, Matrix &R, const Matrix &Z, +InitializeKalmanFilter::initialize(Vector &steadyState, const Vector &deepParams, Matrix &R, const Matrix &Q, Matrix &RQRt, Matrix &T, double &penalty, const MatrixView &dataView, Matrix &Y, int &info) { - modelSolution.compute(steadyState, deepParams, ghx_raw, ghu_raw); + modelSolution.compute(steadyState, deepParams, g_x, g_u); detrendData.detrend(steadyState, dataView, Y); - mat::reorderRowsByVectors(ghx, mat::nullVec, ghx_raw, order_var); - mat::reorderRowsByVectors(ghu, mat::nullVec, ghu_raw, order_var); - setT(T, info); setRQR(R, Q, RQRt, info); } // initialise all KF matrices void -InitializeKalmanFilter::initialize(Vector &steadyState, const Vector &deepParams, Matrix &R, const Matrix &Z, +InitializeKalmanFilter::initialize(Vector &steadyState, const Vector &deepParams, Matrix &R, 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); + initialize(steadyState, deepParams, R, Q, RQRt, T, penalty, dataView, Y, info); setPstar(Pstar, Pinf, T, RQRt, info); } @@ -80,26 +82,19 @@ InitializeKalmanFilter::initialize(Vector &steadyState, const Vector &deepParams void InitializeKalmanFilter::setT(Matrix &T, int &info) { - // here is the content of [T R]=[A,B] = kalman_transition_matrix(oo_.dr,iv,ic,aux,M_.exo_nbr); - + // Initialize the empty columns of T to zero T.setAll(0.0); - - //T(i_n_iv,ic) = dr.ghx(iv,:); - mat::assignByVectors(T, mat::nullVec, inv_ric, ghx, riv, mat::nullVec); + mat::assignByVectors(T, mat::nullVec, pi_bm_vbm, g_x, zeta_varobs_back_mixed, mat::nullVec); } void 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); + mat::assignByVectors(R, mat::nullVec, mat::nullVec, g_u, zeta_varobs_back_mixed, mat::nullVec); // 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' + blas::gemm("N", "N", 1.0, R, Q, 0.0, RQ); // R*Q + blas::gemm("N", "T", 1.0, RQ, R, 0.0, RQRt); // R*Q*R' } void diff --git a/mex/sources/estimation/InitializeKalmanFilter.hh b/mex/sources/estimation/InitializeKalmanFilter.hh index 27af85c8d..982f6a326 100644 --- a/mex/sources/estimation/InitializeKalmanFilter.hh +++ b/mex/sources/estimation/InitializeKalmanFilter.hh @@ -40,29 +40,32 @@ class InitializeKalmanFilter { public: + /*! + \param[in] zeta_varobs_back_mixed_arg The union of indices of observed, backward and mixed variables + */ InitializeKalmanFilter(const std::string& modName, size_t n_endo, size_t n_exo, const std::vector &zeta_fwrd_arg, const std::vector &zeta_back_arg, const std::vector &zeta_mixed_arg, const std::vector &zeta_static_arg, - double qz_criterium, const std::vector &order_var_arg, const std::vector &inv_order_var_arg, - const std::vector &riv, const std::vector &ric, double lyapunov_tol, int &info); + const std::vector &zeta_varobs_back_mixed_arg, + double qz_criterium_arg, double lyapunov_tol_arg, 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, + void initialize(Vector &steadyState, const Vector &deepParams, Matrix &R, 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, + void initialize(Vector &steadyState, const Vector &deepParams, Matrix &R, const Matrix &Q, Matrix &RQRt, Matrix &T, double &penalty, const MatrixView &dataView, Matrix &Y, int &info); private: - const std::vector riv; // restrict_var_list - std::vector inv_ric; // inverse restrict_columns - const std::vector order_var; const double lyapunov_tol; + const std::vector zeta_varobs_back_mixed; + //! Indices of back+mixed zetas inside varobs+back+mixed zetas + std::vector pi_bm_vbm; DetrendData detrendData; ModelSolution modelSolution; DiscLyapFast discLyapFast; //Lyapunov solver - Matrix ghx, ghx_raw; - Matrix ghu, ghu_raw; + Matrix g_x; + Matrix g_u; Matrix Rt, RQ; void setT(Matrix &T, int &info); void setRQR(Matrix &R, const Matrix &Q, Matrix &RQRt, int &info); diff --git a/mex/sources/estimation/KalmanFilter.cc b/mex/sources/estimation/KalmanFilter.cc index 14639da4f..679993e98 100644 --- a/mex/sources/estimation/KalmanFilter.cc +++ b/mex/sources/estimation/KalmanFilter.cc @@ -33,24 +33,41 @@ KalmanFilter::~KalmanFilter() KalmanFilter::KalmanFilter(const std::string &modName, size_t n_endo, size_t n_exo, const std::vector &zeta_fwrd_arg, const std::vector &zeta_back_arg, const std::vector &zeta_mixed_arg, const std::vector &zeta_static_arg, - double qz_criterium, const std::vector &order_var, - const std::vector &inv_order_var, const std::vector &varobs_arg, - const std::vector &riv, const std::vector &ric, - double riccati_tol_in, double lyapunov_tol, int &info) : - Z(varobs_arg.size(), riv.size()), T(riv.size()), R(riv.size(), n_exo), - Pstar(riv.size(), riv.size()), Pinf(riv.size(), riv.size()), RQRt(riv.size(), riv.size()), - Ptmp(riv.size(), riv.size()), F(varobs_arg.size(), varobs_arg.size()), - Finv(varobs_arg.size(), varobs_arg.size()), K(riv.size(), varobs_arg.size()), KFinv(riv.size(), varobs_arg.size()), - oldKFinv(riv.size(), varobs_arg.size()), Finverter(varobs_arg.size()), a_init(riv.size(), 1), - a_new(riv.size(), 1), vt(varobs_arg.size(), 1), vtFinv(1, varobs_arg.size()), vtFinvVt(1), riccati_tol(riccati_tol_in), + double qz_criterium_arg, const std::vector &varobs_arg, + double riccati_tol_arg, double lyapunov_tol_arg, int &info) : + zeta_varobs_back_mixed(compute_zeta_varobs_back_mixed(zeta_back_arg, zeta_mixed_arg, varobs_arg)), + Z(varobs_arg.size(), zeta_varobs_back_mixed.size()), T(zeta_varobs_back_mixed.size()), R(zeta_varobs_back_mixed.size(), n_exo), + Pstar(zeta_varobs_back_mixed.size(), zeta_varobs_back_mixed.size()), Pinf(zeta_varobs_back_mixed.size(), zeta_varobs_back_mixed.size()), RQRt(zeta_varobs_back_mixed.size(), zeta_varobs_back_mixed.size()), + Ptmp(zeta_varobs_back_mixed.size(), zeta_varobs_back_mixed.size()), F(varobs_arg.size(), varobs_arg.size()), + Finv(varobs_arg.size(), varobs_arg.size()), K(zeta_varobs_back_mixed.size(), varobs_arg.size()), KFinv(zeta_varobs_back_mixed.size(), varobs_arg.size()), + oldKFinv(zeta_varobs_back_mixed.size(), varobs_arg.size()), Finverter(varobs_arg.size()), a_init(zeta_varobs_back_mixed.size(), 1), + a_new(zeta_varobs_back_mixed.size(), 1), vt(varobs_arg.size(), 1), vtFinv(1, varobs_arg.size()), vtFinvVt(1), riccati_tol(riccati_tol_arg), initKalmanFilter(modName, n_endo, n_exo, zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg, - zeta_static_arg, qz_criterium, order_var, inv_order_var, riv, ric, lyapunov_tol, info) + zeta_static_arg, zeta_varobs_back_mixed, qz_criterium_arg, lyapunov_tol_arg, info) { a_init.setAll(0.0); Z.setAll(0.0); for (size_t i = 0; i < varobs_arg.size(); ++i) - Z(i, varobs_arg[i]) = 1.0; + { + size_t j = find(zeta_varobs_back_mixed.begin(), zeta_varobs_back_mixed.end(), + varobs_arg[i]) - zeta_varobs_back_mixed.begin(); + Z(i, j) = 1.0; + } +} +std::vector +KalmanFilter::compute_zeta_varobs_back_mixed(const std::vector &zeta_back_arg, const std::vector &zeta_mixed_arg, const std::vector &varobs_arg) +{ + std::vector varobs_sorted = varobs_arg; + sort(varobs_sorted.begin(), varobs_sorted.end()); + std::vector zeta_back_mixed, zeta_varobs_back_mixed; + set_union(zeta_back_arg.begin(), zeta_back_arg.end(), + zeta_mixed_arg.begin(), zeta_mixed_arg.end(), + back_inserter(zeta_back_mixed)); + set_union(zeta_back_mixed.begin(), zeta_back_mixed.end(), + varobs_arg.begin(), varobs_arg.end(), + back_inserter(zeta_varobs_back_mixed)); + return zeta_varobs_back_mixed; } double @@ -58,17 +75,16 @@ KalmanFilter::compute(const MatrixView &dataView, Vector &steadyState, const Matrix &Q, const Matrix &H, const Vector &deepParams, VectorView &vll, size_t start, size_t period, double &penalty, int &info) { - double ll; Matrix Y(dataView.getRows(), dataView.getCols()); // data if(period==0) // initialise all KF matrices - initKalmanFilter.initialize(steadyState, deepParams, R, Z, Q, RQRt, T, Pstar, Pinf, + initKalmanFilter.initialize(steadyState, deepParams, R, 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, + initKalmanFilter.initialize(steadyState, deepParams, R, Q, RQRt, T, penalty, dataView, Y, info); - return ll = filter(Y, H, vll, start, info); + return filter(Y, H, vll, start, info); }; @@ -122,7 +138,6 @@ KalmanFilter::filter(const Matrix &dataView, const Matrix &H, VectorView &vll, } // err= Yt - Za - //VectorView yt(dataView.getData()+t*dataView.getRows(),dataView.getRows(),1); // current observation vector MatrixConstView yt(dataView, 0, t, dataView.getRows(), 1); // current observation vector vt = yt; blas::gemm("N", "N", -1.0, Z, a_init, 1.0, vt); diff --git a/mex/sources/estimation/KalmanFilter.hh b/mex/sources/estimation/KalmanFilter.hh index 4edaa193d..3ccde1068 100644 --- a/mex/sources/estimation/KalmanFilter.hh +++ b/mex/sources/estimation/KalmanFilter.hh @@ -29,7 +29,7 @@ #include "InitializeKalmanFilter.hh" /** - * •vanilla Kalman filter without constant and with measurement error (use scalar + * Vanilla Kalman filter without constant and with measurement error (use scalar * 0 when no measurement error). * If multivariate filter is faster, do as in Matlab: start with multivariate * filter and switch to univariate filter only in case of singularity @@ -50,18 +50,16 @@ public: virtual ~KalmanFilter(); KalmanFilter(const std::string &modName, size_t n_endo, size_t n_exo, const std::vector &zeta_fwrd_arg, const std::vector &zeta_back_arg, const std::vector &zeta_mixed_arg, const std::vector &zeta_static_arg, - double qz_criterium, const std::vector &order_var_arg, const std::vector &inv_order_var, - const std::vector &varobs_arg, const std::vector &riv, const std::vector &ric, - double riccati_tol, double lyapunov_tol, int &info); + double qz_criterium_arg, const std::vector &varobs_arg, + double riccati_tol_arg, double lyapunov_tol_arg, int &info); double compute(const MatrixView &dataView, Vector &steadyState, const Matrix &Q, const Matrix &H, const Vector &deepParams, VectorView &vll, size_t start, size_t period, double &penalty, int &info); -protected: -// static double calcStepLogLik(const PLUFact &Finv, const Vector &v); - private: + const std::vector zeta_varobs_back_mixed; + static std::vector compute_zeta_varobs_back_mixed(const std::vector &zeta_back_arg, const std::vector &zeta_mixed_arg, const std::vector &varobs_arg); Matrix Z; //nob*mm matrix mapping endogeneous variables and observations Matrix T; //mm*mm transition matrix of the state equation. Matrix R; //mm*rr matrix, mapping structural innovations to state variables. diff --git a/mex/sources/estimation/LogLikelihoodMain.hh b/mex/sources/estimation/LogLikelihoodMain.hh index 4f6582333..411bff6f2 100644 --- a/mex/sources/estimation/LogLikelihoodMain.hh +++ b/mex/sources/estimation/LogLikelihoodMain.hh @@ -47,9 +47,9 @@ public: LogLikelihoodMain(const Matrix &data, //GeneralParams& estimOptions, const std::string &modName, EstimatedParametersDescription &estiParDesc, size_t n_endo, size_t n_exo, const std::vector &zeta_fwrd_arg, const std::vector &zeta_back_arg, const std::vector &zeta_mixed_arg, - const std::vector &zeta_static_arg, const double qz_criterium, const std::vector &order_var, - const std::vector &inv_order_var, const std::vector &varobs, const std::vector &riv, - const std::vector &ric, double riccati_tol_in, double lyapunov_tol, int &info); + const std::vector &zeta_static_arg, const double qz_criterium_arg, + const std::vector &varobs_arg, + double riccati_tol_arg, double lyapunov_tol_arg, int &info); double compute(Matrix &steadyState, Vector &estParams, int &info); // for calls from estimation and to set Steady State double compute(Vector &estParams); // for calls in loop from optimiser diff --git a/mex/sources/estimation/LogLikelihoodSubSample.cc b/mex/sources/estimation/LogLikelihoodSubSample.cc index 4193aae9d..7e05c7c56 100644 --- a/mex/sources/estimation/LogLikelihoodSubSample.cc +++ b/mex/sources/estimation/LogLikelihoodSubSample.cc @@ -31,11 +31,9 @@ LogLikelihoodSubSample::~LogLikelihoodSubSample() LogLikelihoodSubSample::LogLikelihoodSubSample(const std::string &modName, size_t n_endo, size_t n_exo, const std::vector &zeta_fwrd_arg, const std::vector &zeta_back_arg, const std::vector &zeta_mixed_arg, - const std::vector &zeta_static_arg, const double qz_criterium, const std::vector &order_var, - const std::vector &inv_order_var, const std::vector &varobs, const std::vector &riv, - const std::vector &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, qz_criterium, - order_var, inv_order_var, varobs, riv, ric, riccati_tol, lyapunov_tol, info) + const std::vector &zeta_static_arg, const double qz_criterium_arg, const std::vector &varobs_arg, double riccati_tol_arg, double lyapunov_tol_arg, int &info) : + kalmanFilter(modName, 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) { }; diff --git a/mex/sources/estimation/LogLikelihoodSubSample.hh b/mex/sources/estimation/LogLikelihoodSubSample.hh index 23b9146c2..0f847e4ac 100644 --- a/mex/sources/estimation/LogLikelihoodSubSample.hh +++ b/mex/sources/estimation/LogLikelihoodSubSample.hh @@ -35,9 +35,8 @@ class LogLikelihoodSubSample { public: LogLikelihoodSubSample(const std::string &modName, EstimatedParametersDescription &estiParDesc, size_t n_endo, size_t n_exo, const std::vector &zeta_fwrd_arg, const std::vector &zeta_back_arg, const std::vector &zeta_mixed_arg, - const std::vector &zeta_static_arg, const double qz_criterium, const std::vector &order_var, - const std::vector &inv_order_var, const std::vector &varobs, const std::vector &riv, - const std::vector &ric, double riccati_tol_in, double lyapunov_tol, int &info); + const std::vector &zeta_static_arg, const double qz_criterium, + const std::vector &varobs_arg, double riccati_tol_in, double lyapunov_tol, int &info); double compute(Vector &steadyState, const MatrixView &dataView, Vector &deepParams, Matrix &Q, Matrix &H, VectorView &vll, int &info, size_t start, size_t period); diff --git a/mex/sources/estimation/tests/testInitKalman.cc b/mex/sources/estimation/tests/testInitKalman.cc index ef7981f6d..c261f9b71 100644 --- a/mex/sources/estimation/tests/testInitKalman.cc +++ b/mex/sources/estimation/tests/testInitKalman.cc @@ -100,47 +100,32 @@ main(int argc, char **argv) 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}; + size_t varobs[] = {12, 11}; + std::vector varobs_arg; + for (size_t i = 0; i < 2; ++i) + varobs_arg.push_back(varobs[i]-1); - std::vector riv; - for (size_t i = 0; i < nriv; ++i) - riv.push_back(sriv[i]-1); - std::vector ric; - for (size_t i = 0; i < nric; ++i) - ric.push_back(sric[i]-1); + // Compute zeta_varobs_back_mixed + sort(varobs_arg.begin(), varobs_arg.end()); + std::vector zeta_back_mixed, zeta_varobs_back_mixed; + set_union(zeta_back_arg.begin(), zeta_back_arg.end(), + zeta_mixed_arg.begin(), zeta_mixed_arg.end(), + back_inserter(zeta_back_mixed)); + set_union(zeta_back_mixed.begin(), zeta_back_mixed.end(), + varobs_arg.begin(), varobs_arg.end(), + back_inserter(zeta_varobs_back_mixed)); - size_t nobs = 2; - size_t svarobs[] = {2, 1}; - std::vector varobs; - for (size_t i = 0; i < nobs; ++i) - varobs.push_back(svarobs[i]-1); + size_t n_vbm = zeta_varobs_back_mixed.size(); Matrix - 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()), Q(n_exo); - Z.setAll(0.0); - for (size_t i = 0; i < varobs.size(); ++i) - Z(i, varobs[i]) = 1.0; + T(n_vbm, n_vbm), R(n_vbm, n_exo), + RQRt(n_vbm, n_vbm), Pstar(n_vbm, n_vbm), + Pinf(n_vbm, n_vbm), Q(n_exo); MatrixView vCovVW(vcov, n_exo, n_exo, n_exo); Q = vCovVW; - size_t sorder_var[] = - { 4, 5, 6, 8, 9, 10, 11, 12, 14, 1, 7, 13, 2, 3, 15}; - std::vector 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 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; int info = 0; const MatrixView @@ -153,13 +138,13 @@ main(int argc, char **argv) InitializeKalmanFilter initializeKalmanFilter(modName, n_endo, n_exo, - zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg, zeta_static_arg, qz_criterium, - order_var, inv_order_var, riv, ric, lyapunov_tol, info); + zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg, zeta_static_arg, + zeta_varobs_back_mixed, qz_criterium, + lyapunov_tol, info); - std::cout << "Initilise KF with Q: " << std::endl << Q << std::endl; - std::cout << "and Z" << std::endl << Z << std::endl; + std::cout << "Initialize KF with Q: " << std::endl << Q << std::endl; - initializeKalmanFilter.initialize(steadyState, deepParams, R, Z, Q, RQRt, T, Pstar, Pinf, + initializeKalmanFilter.initialize(steadyState, deepParams, R, Q, RQRt, T, Pstar, Pinf, penalty, dataView, yView, info); std::cout << "Matrix T: " << std::endl << T << std::endl; diff --git a/mex/sources/estimation/tests/testKalman.cc b/mex/sources/estimation/tests/testKalman.cc index 947690fdb..b240a1bce 100644 --- a/mex/sources/estimation/tests/testKalman.cc +++ b/mex/sources/estimation/tests/testKalman.cc @@ -100,51 +100,20 @@ main(int argc, char **argv) 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 riv; - for (size_t i = 0; i < nriv; ++i) - riv.push_back(sriv[i]-1); - std::vector 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 varobs; + size_t varobs[] = {12, 11}; + std::vector varobs_arg; for (size_t i = 0; i < nobs; ++i) - varobs.push_back(svarobs[i]-1); + varobs_arg.push_back(varobs[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 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 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; @@ -158,10 +127,8 @@ main(int argc, char **argv) KalmanFilter kalman(modName, n_endo, n_exo, zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg, zeta_static_arg, qz_criterium, - order_var, inv_order_var, varobs, riv, ric, riccati_tol, lyapunov_tol, info); + varobs_arg, 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, period=0; double ll=kalman.compute(dataView, steadyState, Q, H, deepParams, vwll, start, period, penalty, info);