Estimation DLL: simplification of the handling of indices in T and R matrices

time-shift
Sébastien Villemot 2010-05-06 18:57:15 +02:00
parent 5aa17370e8
commit e4f640f201
9 changed files with 112 additions and 152 deletions

View File

@ -33,46 +33,48 @@ InitializeKalmanFilter::~InitializeKalmanFilter()
InitializeKalmanFilter::InitializeKalmanFilter(const std::string &modName, size_t n_endo_arg, size_t n_exo_arg, InitializeKalmanFilter::InitializeKalmanFilter(const std::string &modName, size_t n_endo_arg, size_t n_exo_arg,
const std::vector<size_t> &zeta_fwrd_arg, const std::vector<size_t> &zeta_back_arg, 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 std::vector<size_t> &zeta_mixed_arg, const std::vector<size_t> &zeta_static_arg,
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> &zeta_varobs_back_mixed_arg,
const std::vector<size_t> &rivIN, const std::vector<size_t> &ricIN, double lyapunov_tolIN, int &info) : double qz_criterium_arg,
riv(rivIN), inv_ric(ricIN.size()), order_var(order_var_arg), lyapunov_tol(lyapunov_tolIN), 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, 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()), zeta_mixed_arg, zeta_static_arg, qz_criterium_arg),
ghx(n_endo_arg, zeta_back_arg.size() + zeta_mixed_arg.size()), discLyapFast(zeta_varobs_back_mixed.size()),
ghx_raw(n_endo_arg, zeta_back_arg.size() + zeta_mixed_arg.size()), g_x(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), g_u(n_endo_arg, n_exo_arg),
Rt(n_exo_arg, riv.size()), RQ(riv.size(), 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(); std::vector<size_t> zeta_back_mixed;
size_t n_static = zeta_static_arg.size(); set_union(zeta_back_arg.begin(), zeta_back_arg.end(),
size_t j = 0; zeta_mixed_arg.begin(), zeta_mixed_arg.end(),
for (size_t i = 0; i < n_endo_arg; ++i) back_inserter(zeta_back_mixed));
if (inv_order_var_arg[i]-n_static < n_pred && inv_order_var_arg[i]-n_static >= 0) for (size_t i = 0; i < zeta_back_mixed.size(); i++)
inv_ric[j++] = ricIN[inv_order_var_arg[i]-n_static]; 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 // initialise parameter dependent KF matrices only but not Ps
void 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, const Matrix &Q, Matrix &RQRt, Matrix &T,
double &penalty, const MatrixView &dataView, Matrix &Y, int &info) 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); 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); setT(T, info);
setRQR(R, Q, RQRt, info); setRQR(R, Q, RQRt, info);
} }
// initialise all KF matrices // initialise all KF matrices
void 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, const Matrix &Q, Matrix &RQRt, Matrix &T, Matrix &Pstar, Matrix &Pinf,
double &penalty, const MatrixView &dataView, Matrix &Y, int &info) 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); setPstar(Pstar, Pinf, T, RQRt, info);
} }
@ -80,26 +82,19 @@ InitializeKalmanFilter::initialize(Vector &steadyState, const Vector &deepParams
void void
InitializeKalmanFilter::setT(Matrix &T, int &info) 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.setAll(0.0);
mat::assignByVectors(T, mat::nullVec, pi_bm_vbm, g_x, zeta_varobs_back_mixed, mat::nullVec);
//T(i_n_iv,ic) = dr.ghx(iv,:);
mat::assignByVectors(T, mat::nullVec, inv_ric, ghx, riv, mat::nullVec);
} }
void void
InitializeKalmanFilter::setRQR(Matrix &R, const Matrix &Q, Matrix &RQRt, int &info) InitializeKalmanFilter::setRQR(Matrix &R, const Matrix &Q, Matrix &RQRt, int &info)
{ {
R.setAll(0.0); mat::assignByVectors(R, mat::nullVec, mat::nullVec, g_u, zeta_varobs_back_mixed, mat::nullVec);
//B(i_n_iv,:) = dr.ghu(iv,:); and R=B;
mat::assignByVectors(R, mat::nullVec, mat::nullVec, ghu, riv, mat::nullVec);
// Matrix RQRt=R*Q*R' // 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 blas::gemm("N", "T", 1.0, RQ, R, 0.0, RQRt); // R*Q*R'
RQRt.setAll(0.0);
blas::gemm("N", "T", 1.0, RQ, R, 1.0, RQRt); // R*Q*R'
} }
void void

View File

@ -40,29 +40,32 @@ class InitializeKalmanFilter
{ {
public: 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<size_t> &zeta_fwrd_arg, InitializeKalmanFilter(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 std::vector<size_t> &zeta_back_arg, const std::vector<size_t> &zeta_mixed_arg, const std::vector<size_t> &zeta_static_arg,
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> &zeta_varobs_back_mixed_arg,
const std::vector<size_t> &riv, const std::vector<size_t> &ric, double lyapunov_tol, int &info); double qz_criterium_arg, double lyapunov_tol_arg, int &info);
virtual ~InitializeKalmanFilter(); virtual ~InitializeKalmanFilter();
// initialise all KF matrices // 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); Matrix &T, Matrix &Pstar, Matrix &Pinf, double &penalty, const MatrixView &dataView, Matrix &Y, int &info);
// initialise parameter dependent KF matrices only but not Ps // 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); Matrix &T, double &penalty, const MatrixView &dataView, Matrix &Y, int &info);
private: private:
const std::vector<size_t> riv; // restrict_var_list
std::vector<size_t> inv_ric; // inverse restrict_columns
const std::vector<size_t> order_var;
const double lyapunov_tol; const double lyapunov_tol;
const std::vector<size_t> zeta_varobs_back_mixed;
//! Indices of back+mixed zetas inside varobs+back+mixed zetas
std::vector<size_t> pi_bm_vbm;
DetrendData detrendData; DetrendData detrendData;
ModelSolution modelSolution; ModelSolution modelSolution;
DiscLyapFast discLyapFast; //Lyapunov solver DiscLyapFast discLyapFast; //Lyapunov solver
Matrix ghx, ghx_raw; Matrix g_x;
Matrix ghu, ghu_raw; Matrix g_u;
Matrix Rt, RQ; Matrix Rt, RQ;
void setT(Matrix &T, int &info); void setT(Matrix &T, int &info);
void setRQR(Matrix &R, const Matrix &Q, Matrix &RQRt, int &info); void setRQR(Matrix &R, const Matrix &Q, Matrix &RQRt, int &info);

View File

@ -33,24 +33,41 @@ KalmanFilter::~KalmanFilter()
KalmanFilter::KalmanFilter(const std::string &modName, size_t n_endo, size_t n_exo, KalmanFilter::KalmanFilter(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_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 std::vector<size_t> &zeta_mixed_arg, const std::vector<size_t> &zeta_static_arg,
double qz_criterium, const std::vector<size_t> &order_var, double qz_criterium_arg, const std::vector<size_t> &varobs_arg,
const std::vector<size_t> &inv_order_var, const std::vector<size_t> &varobs_arg, double riccati_tol_arg, double lyapunov_tol_arg, int &info) :
const std::vector<size_t> &riv, const std::vector<size_t> &ric, zeta_varobs_back_mixed(compute_zeta_varobs_back_mixed(zeta_back_arg, zeta_mixed_arg, varobs_arg)),
double riccati_tol_in, double lyapunov_tol, int &info) : Z(varobs_arg.size(), zeta_varobs_back_mixed.size()), T(zeta_varobs_back_mixed.size()), R(zeta_varobs_back_mixed.size(), n_exo),
Z(varobs_arg.size(), riv.size()), T(riv.size()), R(riv.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()),
Pstar(riv.size(), riv.size()), Pinf(riv.size(), riv.size()), RQRt(riv.size(), riv.size()), Ptmp(zeta_varobs_back_mixed.size(), zeta_varobs_back_mixed.size()), F(varobs_arg.size(), varobs_arg.size()),
Ptmp(riv.size(), riv.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()),
Finv(varobs_arg.size(), varobs_arg.size()), K(riv.size(), varobs_arg.size()), KFinv(riv.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),
oldKFinv(riv.size(), varobs_arg.size()), Finverter(varobs_arg.size()), a_init(riv.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),
a_new(riv.size(), 1), vt(varobs_arg.size(), 1), vtFinv(1, varobs_arg.size()), vtFinvVt(1), riccati_tol(riccati_tol_in),
initKalmanFilter(modName, n_endo, n_exo, zeta_fwrd_arg, zeta_back_arg, zeta_mixed_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); a_init.setAll(0.0);
Z.setAll(0.0); Z.setAll(0.0);
for (size_t i = 0; i < varobs_arg.size(); ++i) 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<size_t>
KalmanFilter::compute_zeta_varobs_back_mixed(const std::vector<size_t> &zeta_back_arg, const std::vector<size_t> &zeta_mixed_arg, const std::vector<size_t> &varobs_arg)
{
std::vector<size_t> varobs_sorted = varobs_arg;
sort(varobs_sorted.begin(), varobs_sorted.end());
std::vector<size_t> 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 double
@ -58,17 +75,16 @@ KalmanFilter::compute(const MatrixView &dataView, Vector &steadyState,
const Matrix &Q, const Matrix &H, const Vector &deepParams, const Matrix &Q, const Matrix &H, const Vector &deepParams,
VectorView &vll, size_t start, size_t period, 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 Matrix Y(dataView.getRows(), dataView.getCols()); // data
if(period==0) // initialise all KF matrices 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); penalty, dataView, Y, info);
else // initialise parameter dependent KF matrices only but not Ps 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); 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 // 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 MatrixConstView yt(dataView, 0, t, dataView.getRows(), 1); // current observation vector
vt = yt; vt = yt;
blas::gemm("N", "N", -1.0, Z, a_init, 1.0, vt); blas::gemm("N", "N", -1.0, Z, a_init, 1.0, vt);

View File

@ -29,7 +29,7 @@
#include "InitializeKalmanFilter.hh" #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). * 0 when no measurement error).
* If multivariate filter is faster, do as in Matlab: start with multivariate * If multivariate filter is faster, do as in Matlab: start with multivariate
* filter and switch to univariate filter only in case of singularity * filter and switch to univariate filter only in case of singularity
@ -50,18 +50,16 @@ public:
virtual ~KalmanFilter(); virtual ~KalmanFilter();
KalmanFilter(const std::string &modName, size_t n_endo, size_t n_exo, const std::vector<size_t> &zeta_fwrd_arg, KalmanFilter(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 std::vector<size_t> &zeta_back_arg, const std::vector<size_t> &zeta_mixed_arg, const std::vector<size_t> &zeta_static_arg,
double qz_criterium, const std::vector<size_t> &order_var_arg, const std::vector<size_t> &inv_order_var, double qz_criterium_arg, const std::vector<size_t> &varobs_arg,
const std::vector<size_t> &varobs_arg, const std::vector<size_t> &riv, const std::vector<size_t> &ric, double riccati_tol_arg, double lyapunov_tol_arg, int &info);
double riccati_tol, double lyapunov_tol, int &info);
double compute(const MatrixView &dataView, Vector &steadyState, double compute(const MatrixView &dataView, Vector &steadyState,
const Matrix &Q, const Matrix &H, const Vector &deepParams, const Matrix &Q, const Matrix &H, const Vector &deepParams,
VectorView &vll, size_t start, size_t period, 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);
private: private:
const std::vector<size_t> zeta_varobs_back_mixed;
static std::vector<size_t> compute_zeta_varobs_back_mixed(const std::vector<size_t> &zeta_back_arg, const std::vector<size_t> &zeta_mixed_arg, const std::vector<size_t> &varobs_arg);
Matrix Z; //nob*mm matrix mapping endogeneous variables and observations Matrix Z; //nob*mm matrix mapping endogeneous variables and observations
Matrix T; //mm*mm transition matrix of the state equation. Matrix T; //mm*mm transition matrix of the state equation.
Matrix R; //mm*rr matrix, mapping structural innovations to state variables. Matrix R; //mm*rr matrix, mapping structural innovations to state variables.

View File

@ -47,9 +47,9 @@ public:
LogLikelihoodMain(const Matrix &data, //GeneralParams& estimOptions, LogLikelihoodMain(const Matrix &data, //GeneralParams& estimOptions,
const std::string &modName, EstimatedParametersDescription &estiParDesc, size_t n_endo, size_t n_exo, const std::string &modName, EstimatedParametersDescription &estiParDesc, 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_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> &order_var, const std::vector<size_t> &zeta_static_arg, const double qz_criterium_arg,
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> &varobs_arg,
const std::vector<size_t> &ric, double riccati_tol_in, double lyapunov_tol, int &info); 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(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 double compute(Vector &estParams); // for calls in loop from optimiser

View File

@ -31,11 +31,9 @@ LogLikelihoodSubSample::~LogLikelihoodSubSample()
LogLikelihoodSubSample::LogLikelihoodSubSample(const std::string &modName, size_t n_endo, size_t n_exo, 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_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> &order_var, const std::vector<size_t> &zeta_static_arg, const double qz_criterium_arg, const std::vector<size_t> &varobs_arg, double riccati_tol_arg, double lyapunov_tol_arg, int &info) :
const std::vector<size_t> &inv_order_var, const std::vector<size_t> &varobs, const std::vector<size_t> &riv, kalmanFilter(modName, n_endo, n_exo, zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg, zeta_static_arg, qz_criterium_arg,
const std::vector<size_t> &ric, double riccati_tol, double lyapunov_tol, int &info) : varobs_arg, riccati_tol_arg, lyapunov_tol_arg, 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)
{ {
}; };

View File

@ -35,9 +35,8 @@ class LogLikelihoodSubSample {
public: public:
LogLikelihoodSubSample(const std::string &modName, EstimatedParametersDescription &estiParDesc, size_t n_endo, size_t n_exo, LogLikelihoodSubSample(const std::string &modName, EstimatedParametersDescription &estiParDesc, 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_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> &order_var, const std::vector<size_t> &zeta_static_arg, const double qz_criterium,
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> &varobs_arg, double riccati_tol_in, double lyapunov_tol, int &info);
const std::vector<size_t> &ric, double riccati_tol_in, double lyapunov_tol, int &info);
double compute(Vector &steadyState, const MatrixView &dataView, Vector &deepParams, double compute(Vector &steadyState, const MatrixView &dataView, Vector &deepParams,
Matrix &Q, Matrix &H, VectorView &vll, int &info, size_t start, size_t period); Matrix &Q, Matrix &H, VectorView &vll, int &info, size_t start, size_t period);

View File

@ -100,47 +100,32 @@ main(int argc, char **argv)
for (int i = 0; i < 2; ++i) for (int i = 0; i < 2; ++i)
zeta_fwrd_arg.push_back(fwd[i]-1); zeta_fwrd_arg.push_back(fwd[i]-1);
size_t nriv = 6, nric = 4; size_t varobs[] = {12, 11};
size_t sriv[] = {7, 8, 10, 11, 12, 13}; std::vector<size_t> varobs_arg;
size_t sric[] = {3, 4, 5, 6}; for (size_t i = 0; i < 2; ++i)
varobs_arg.push_back(varobs[i]-1);
std::vector<size_t> riv; // Compute zeta_varobs_back_mixed
for (size_t i = 0; i < nriv; ++i) sort(varobs_arg.begin(), varobs_arg.end());
riv.push_back(sriv[i]-1); std::vector<size_t> zeta_back_mixed, zeta_varobs_back_mixed;
std::vector<size_t> ric; set_union(zeta_back_arg.begin(), zeta_back_arg.end(),
for (size_t i = 0; i < nric; ++i) zeta_mixed_arg.begin(), zeta_mixed_arg.end(),
ric.push_back(sric[i]-1); 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 n_vbm = zeta_varobs_back_mixed.size();
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 Matrix
T(riv.size(), riv.size()), R(riv.size(), n_exo), T(n_vbm, n_vbm), R(n_vbm, n_exo),
RQRt(riv.size(), riv.size()), Pstar(riv.size(), riv.size()), RQRt(n_vbm, n_vbm), Pstar(n_vbm, n_vbm),
Pinf(riv.size(), riv.size()), Z(varobs.size(), riv.size()), Q(n_exo); Pinf(n_vbm, n_vbm), Q(n_exo);
Z.setAll(0.0);
for (size_t i = 0; i < varobs.size(); ++i)
Z(i, varobs[i]) = 1.0;
MatrixView MatrixView
vCovVW(vcov, n_exo, n_exo, n_exo); vCovVW(vcov, n_exo, n_exo, n_exo);
Q = vCovVW; Q = vCovVW;
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 lyapunov_tol = 1e-16;
int info = 0; int info = 0;
const MatrixView const MatrixView
@ -153,13 +138,13 @@ main(int argc, char **argv)
InitializeKalmanFilter InitializeKalmanFilter
initializeKalmanFilter(modName, n_endo, n_exo, initializeKalmanFilter(modName, n_endo, n_exo,
zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg, zeta_static_arg, qz_criterium, zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg, zeta_static_arg,
order_var, inv_order_var, riv, ric, lyapunov_tol, info); zeta_varobs_back_mixed, qz_criterium,
lyapunov_tol, info);
std::cout << "Initilise KF with Q: " << std::endl << Q << std::endl; std::cout << "Initialize KF with Q: " << std::endl << Q << std::endl;
std::cout << "and Z" << std::endl << Z << 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); penalty, dataView, yView, info);
std::cout << "Matrix T: " << std::endl << T << std::endl; std::cout << "Matrix T: " << std::endl << T << std::endl;

View File

@ -100,51 +100,20 @@ main(int argc, char **argv)
for (int i = 0; i < 2; ++i) for (int i = 0; i < 2; ++i)
zeta_fwrd_arg.push_back(fwd[i]-1); 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 nobs = 2;
size_t svarobs[] = {2, 1}; size_t varobs[] = {12, 11};
std::vector<size_t> varobs; std::vector<size_t> varobs_arg;
for (size_t i = 0; i < nobs; ++i) 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); 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); H.setAll(0.0);
/* for (size_t i = 0; i < varobs.size(); ++i)
Z(i, varobs[i]) = 1.0;
*/
MatrixView MatrixView
vCovVW(vcov, n_exo, n_exo, n_exo); vCovVW(vcov, n_exo, n_exo, n_exo);
Q = vCovVW; Q = vCovVW;
std::cout << "Matrix Q: " << std::endl << Q << std::endl; 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 lyapunov_tol = 1e-16;
double riccati_tol = 1e-16; double riccati_tol = 1e-16;
int info = 0; int info = 0;
@ -158,10 +127,8 @@ main(int argc, char **argv)
KalmanFilter kalman(modName, n_endo, n_exo, KalmanFilter kalman(modName, n_endo, n_exo,
zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg, zeta_static_arg, qz_criterium, 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; size_t start=0, period=0;
double ll=kalman.compute(dataView, steadyState, Q, H, deepParams, double ll=kalman.compute(dataView, steadyState, Q, H, deepParams,
vwll, start, period, penalty, info); vwll, start, period, penalty, info);