Start refactoring of penalty handling
Handle all events that lead to a penalty through an exception. Next steps: - make all exceptions derive from a common class containing info code and penalty value - catch these exceptions at the top-level, and construct the penalty-adjusted likelihoodtime-shift
parent
82930ee29a
commit
f89edf476e
|
@ -37,8 +37,7 @@ InitializeKalmanFilter::InitializeKalmanFilter(const std::string &basename, size
|
||||||
const std::vector<size_t> &varobs_arg,
|
const std::vector<size_t> &varobs_arg,
|
||||||
double qz_criterium_arg,
|
double qz_criterium_arg,
|
||||||
double lyapunov_tol_arg,
|
double lyapunov_tol_arg,
|
||||||
bool noconstant_arg,
|
bool noconstant_arg) :
|
||||||
int &info) :
|
|
||||||
lyapunov_tol(lyapunov_tol_arg),
|
lyapunov_tol(lyapunov_tol_arg),
|
||||||
zeta_varobs_back_mixed(zeta_varobs_back_mixed_arg),
|
zeta_varobs_back_mixed(zeta_varobs_back_mixed_arg),
|
||||||
detrendData(varobs_arg, noconstant_arg),
|
detrendData(varobs_arg, noconstant_arg),
|
||||||
|
@ -60,7 +59,7 @@ InitializeKalmanFilter::InitializeKalmanFilter(const std::string &basename, size
|
||||||
|
|
||||||
}
|
}
|
||||||
void
|
void
|
||||||
InitializeKalmanFilter::setT(Matrix &T, int &info)
|
InitializeKalmanFilter::setT(Matrix &T)
|
||||||
{
|
{
|
||||||
// Initialize the empty columns of T to zero
|
// Initialize the empty columns of T to zero
|
||||||
T.setAll(0.0);
|
T.setAll(0.0);
|
||||||
|
@ -68,31 +67,11 @@ InitializeKalmanFilter::setT(Matrix &T, int &info)
|
||||||
}
|
}
|
||||||
|
|
||||||
void
|
void
|
||||||
InitializeKalmanFilter::setPstar(Matrix &Pstar, Matrix &Pinf, const Matrix &T, const Matrix &RQRt, int &info)
|
InitializeKalmanFilter::setPstar(Matrix &Pstar, Matrix &Pinf, const Matrix &T, const Matrix &RQRt) throw (DiscLyapFast::DLPException)
|
||||||
{
|
{
|
||||||
|
|
||||||
try
|
|
||||||
{
|
|
||||||
// disclyap_fast(T, RQR, Pstar, lyapunov_tol, 0 or 1 to check chol)
|
|
||||||
discLyapFast.solve_lyap(T, RQRt, Pstar, lyapunov_tol, 0);
|
discLyapFast.solve_lyap(T, RQRt, Pstar, lyapunov_tol, 0);
|
||||||
|
|
||||||
Pinf.setAll(0.0);
|
Pinf.setAll(0.0);
|
||||||
}
|
|
||||||
catch (const DiscLyapFast::DLPException &e)
|
|
||||||
{
|
|
||||||
if (e.info > 0) // The matrix is not positive definite in NormCholesky calculator
|
|
||||||
{
|
|
||||||
puts(e.message.c_str());
|
|
||||||
info = -1; //likelihood = penalty;
|
|
||||||
return;
|
|
||||||
}
|
|
||||||
else if (e.info < 0)
|
|
||||||
{
|
|
||||||
printf("Caugth unhandled TS exception with Pstar matrix: ");
|
|
||||||
puts(e.message.c_str());
|
|
||||||
info = -1; //likelihood = penalty;
|
|
||||||
throw;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
|
@ -48,31 +48,31 @@ public:
|
||||||
const std::vector<size_t> &zeta_varobs_back_mixed_arg,
|
const std::vector<size_t> &zeta_varobs_back_mixed_arg,
|
||||||
const std::vector<size_t> &varobs_arg,
|
const std::vector<size_t> &varobs_arg,
|
||||||
double qz_criterium_arg, double lyapunov_tol_arg,
|
double qz_criterium_arg, double lyapunov_tol_arg,
|
||||||
bool noconstant_arg, int &info);
|
bool noconstant_arg);
|
||||||
virtual ~InitializeKalmanFilter();
|
virtual ~InitializeKalmanFilter();
|
||||||
// initialise parameter dependent KF matrices only but not Ps
|
// initialise parameter dependent KF matrices only but not Ps
|
||||||
template <class Vec1, class Vec2, class Mat1, class Mat2>
|
template <class Vec1, class Vec2, class Mat1, class Mat2>
|
||||||
void initialize(Vec1 &steadyState, const Vec2 &deepParams, Mat1 &R,
|
void initialize(Vec1 &steadyState, const Vec2 &deepParams, Mat1 &R,
|
||||||
const Mat2 &Q, Matrix &RQRt, Matrix &T,
|
const Mat2 &Q, Matrix &RQRt, Matrix &T,
|
||||||
double &penalty, const MatrixConstView &dataView,
|
const MatrixConstView &dataView,
|
||||||
MatrixView &detrendedDataView, int &info)
|
MatrixView &detrendedDataView)
|
||||||
{
|
{
|
||||||
modelSolution.compute(steadyState, deepParams, g_x, g_u);
|
modelSolution.compute(steadyState, deepParams, g_x, g_u);
|
||||||
detrendData.detrend(steadyState, dataView, detrendedDataView);
|
detrendData.detrend(steadyState, dataView, detrendedDataView);
|
||||||
|
|
||||||
setT(T, info);
|
setT(T);
|
||||||
setRQR(R, Q, RQRt, info);
|
setRQR(R, Q, RQRt);
|
||||||
}
|
}
|
||||||
|
|
||||||
// initialise all KF matrices
|
// initialise all KF matrices
|
||||||
template <class Vec1, class Vec2, class Mat1, class Mat2>
|
template <class Vec1, class Vec2, class Mat1, class Mat2>
|
||||||
void initialize(Vec1 &steadyState, const Vec2 &deepParams, Mat1 &R,
|
void initialize(Vec1 &steadyState, const Vec2 &deepParams, Mat1 &R,
|
||||||
const Mat2 &Q, Matrix &RQRt, Matrix &T, Matrix &Pstar, Matrix &Pinf,
|
const Mat2 &Q, Matrix &RQRt, Matrix &T, Matrix &Pstar, Matrix &Pinf,
|
||||||
double &penalty, const MatrixConstView &dataView,
|
const MatrixConstView &dataView,
|
||||||
MatrixView &detrendedDataView, int &info)
|
MatrixView &detrendedDataView)
|
||||||
{
|
{
|
||||||
initialize(steadyState, deepParams, R, Q, RQRt, T, penalty, dataView, detrendedDataView, info);
|
initialize(steadyState, deepParams, R, Q, RQRt, T, dataView, detrendedDataView);
|
||||||
setPstar(Pstar, Pinf, T, RQRt, info);
|
setPstar(Pstar, Pinf, T, RQRt);
|
||||||
}
|
}
|
||||||
|
|
||||||
private:
|
private:
|
||||||
|
@ -87,10 +87,10 @@ private:
|
||||||
Matrix g_x;
|
Matrix g_x;
|
||||||
Matrix g_u;
|
Matrix g_u;
|
||||||
Matrix Rt, RQ;
|
Matrix Rt, RQ;
|
||||||
void setT(Matrix &T, int &info);
|
void setT(Matrix &T);
|
||||||
|
|
||||||
template <class Mat1, class Mat2>
|
template <class Mat1, class Mat2>
|
||||||
void setRQR(Mat1 &R, const Mat2 &Q, Matrix &RQRt, int &info)
|
void setRQR(Mat1 &R, const Mat2 &Q, Matrix &RQRt)
|
||||||
{
|
{
|
||||||
mat::assignByVectors(R, mat::nullVec, mat::nullVec, g_u, zeta_varobs_back_mixed, mat::nullVec);
|
mat::assignByVectors(R, mat::nullVec, mat::nullVec, g_u, zeta_varobs_back_mixed, mat::nullVec);
|
||||||
|
|
||||||
|
@ -98,7 +98,7 @@ private:
|
||||||
blas::gemm("N", "N", 1.0, R, Q, 0.0, RQ); // R*Q
|
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'
|
blas::gemm("N", "T", 1.0, RQ, R, 0.0, RQRt); // R*Q*R'
|
||||||
}
|
}
|
||||||
void setPstar(Matrix &Pstar, Matrix &Pinf, const Matrix &T, const Matrix &RQRt, int &info);
|
void setPstar(Matrix &Pstar, Matrix &Pinf, const Matrix &T, const Matrix &RQRt) throw (DiscLyapFast::DLPException);
|
||||||
|
|
||||||
};
|
};
|
||||||
|
|
||||||
|
|
|
@ -1,5 +1,5 @@
|
||||||
/*
|
/*
|
||||||
* Copyright (C) 2009-2012 Dynare Team
|
* Copyright (C) 2009-2013 Dynare Team
|
||||||
*
|
*
|
||||||
* This file is part of Dynare.
|
* This file is part of Dynare.
|
||||||
*
|
*
|
||||||
|
@ -36,7 +36,7 @@ KalmanFilter::KalmanFilter(const std::string &basename, size_t n_endo, size_t n_
|
||||||
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_arg, const std::vector<size_t> &varobs_arg,
|
double qz_criterium_arg, const std::vector<size_t> &varobs_arg,
|
||||||
double riccati_tol_arg, double lyapunov_tol_arg,
|
double riccati_tol_arg, double lyapunov_tol_arg,
|
||||||
bool noconstant_arg, int &info) :
|
bool noconstant_arg) :
|
||||||
zeta_varobs_back_mixed(compute_zeta_varobs_back_mixed(zeta_back_arg, zeta_mixed_arg, varobs_arg)),
|
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()), Zt(Z.getCols(), Z.getRows()), T(zeta_varobs_back_mixed.size()), R(zeta_varobs_back_mixed.size(), n_exo),
|
Z(varobs_arg.size(), zeta_varobs_back_mixed.size()), Zt(Z.getCols(), Z.getRows()), 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()),
|
Pstar(zeta_varobs_back_mixed.size(), zeta_varobs_back_mixed.size()), Pinf(zeta_varobs_back_mixed.size(), zeta_varobs_back_mixed.size()),
|
||||||
|
@ -45,7 +45,7 @@ KalmanFilter::KalmanFilter(const std::string &basename, size_t n_endo, size_t n_
|
||||||
oldKFinv(zeta_varobs_back_mixed.size(), varobs_arg.size()), a_init(zeta_varobs_back_mixed.size()),
|
oldKFinv(zeta_varobs_back_mixed.size(), varobs_arg.size()), a_init(zeta_varobs_back_mixed.size()),
|
||||||
a_new(zeta_varobs_back_mixed.size()), vt(varobs_arg.size()), vtFinv(varobs_arg.size()), riccati_tol(riccati_tol_arg),
|
a_new(zeta_varobs_back_mixed.size()), vt(varobs_arg.size()), vtFinv(varobs_arg.size()), riccati_tol(riccati_tol_arg),
|
||||||
initKalmanFilter(basename, n_endo, n_exo, zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg,
|
initKalmanFilter(basename, n_endo, n_exo, zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg,
|
||||||
zeta_static_arg, zeta_varobs_back_mixed, varobs_arg, qz_criterium_arg, lyapunov_tol_arg, noconstant_arg, info),
|
zeta_static_arg, zeta_varobs_back_mixed, varobs_arg, qz_criterium_arg, lyapunov_tol_arg, noconstant_arg),
|
||||||
FUTP(varobs_arg.size()*(varobs_arg.size()+1)/2)
|
FUTP(varobs_arg.size()*(varobs_arg.size()+1)/2)
|
||||||
{
|
{
|
||||||
Z.setAll(0.0);
|
Z.setAll(0.0);
|
||||||
|
@ -79,12 +79,14 @@ KalmanFilter::compute_zeta_varobs_back_mixed(const std::vector<size_t> &zeta_bac
|
||||||
* Multi-variate standard Kalman Filter
|
* Multi-variate standard Kalman Filter
|
||||||
*/
|
*/
|
||||||
double
|
double
|
||||||
KalmanFilter::filter(const MatrixView &detrendedDataView, const Matrix &H, VectorView &vll, size_t start, int &info)
|
KalmanFilter::filter(const MatrixView &detrendedDataView, const Matrix &H, VectorView &vll, size_t start)
|
||||||
{
|
{
|
||||||
double loglik = 0.0, ll, logFdet = 0.0, Fdet, dvtFinvVt;
|
double loglik = 0.0, ll, logFdet = 0.0, Fdet, dvtFinvVt;
|
||||||
size_t p = Finv.getRows();
|
size_t p = Finv.getRows();
|
||||||
bool nonstationary = true;
|
bool nonstationary = true;
|
||||||
a_init.setAll(0.0);
|
a_init.setAll(0.0);
|
||||||
|
int info;
|
||||||
|
|
||||||
for (size_t t = 0; t < detrendedDataView.getCols(); ++t)
|
for (size_t t = 0; t < detrendedDataView.getCols(); ++t)
|
||||||
{
|
{
|
||||||
if (nonstationary)
|
if (nonstationary)
|
||||||
|
@ -106,12 +108,8 @@ KalmanFilter::filter(const MatrixView &detrendedDataView, const Matrix &H, Vect
|
||||||
FUTP(i + (j-1)*j/2 -1) = F(i-1, j-1);
|
FUTP(i + (j-1)*j/2 -1) = F(i-1, j-1);
|
||||||
|
|
||||||
info = lapack::choleskySolver(FUTP, Finv, "U"); // F now contains its Chol decomposition!
|
info = lapack::choleskySolver(FUTP, Finv, "U"); // F now contains its Chol decomposition!
|
||||||
if (info < 0)
|
assert(info >= 0);
|
||||||
{
|
|
||||||
std::cout << "Info:" << info << std::endl;
|
|
||||||
std::cout << "t:" << t << std::endl;
|
|
||||||
return 0;
|
|
||||||
}
|
|
||||||
if (info > 0)
|
if (info > 0)
|
||||||
{
|
{
|
||||||
//enforce Pstar symmetry with P=(P+P')/2=0.5P+0.5P'
|
//enforce Pstar symmetry with P=(P+P')/2=0.5P+0.5P'
|
||||||
|
@ -136,11 +134,10 @@ KalmanFilter::filter(const MatrixView &detrendedDataView, const Matrix &H, Vect
|
||||||
for (size_t j = i; j <= p; ++j)
|
for (size_t j = i; j <= p; ++j)
|
||||||
FUTP(i + (j-1)*j/2 -1) = F(i-1, j-1);
|
FUTP(i + (j-1)*j/2 -1) = F(i-1, j-1);
|
||||||
|
|
||||||
info = lapack::choleskySolver(FUTP, Finv, "U"); // F now contains its Chol decomposition!
|
info = lapack::choleskySolver(FUTP, Finv, "U"); // F now contains
|
||||||
if (info != 0)
|
// its Chol
|
||||||
{
|
// decomposition!
|
||||||
return 0;
|
assert(info == 0);
|
||||||
}
|
|
||||||
}
|
}
|
||||||
// KFinv gain matrix
|
// KFinv gain matrix
|
||||||
blas::symm("R", "U", 1.0, Finv, K, 0.0, KFinv);
|
blas::symm("R", "U", 1.0, Finv, K, 0.0, KFinv);
|
||||||
|
|
|
@ -53,38 +53,24 @@ public:
|
||||||
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_arg, const std::vector<size_t> &varobs_arg,
|
double qz_criterium_arg, const std::vector<size_t> &varobs_arg,
|
||||||
double riccati_tol_arg, double lyapunov_tol_arg,
|
double riccati_tol_arg, double lyapunov_tol_arg,
|
||||||
bool noconstant_arg, int &info);
|
bool noconstant_arg);
|
||||||
|
|
||||||
template <class Vec1, class Vec2, class Mat1>
|
template <class Vec1, class Vec2, class Mat1>
|
||||||
double compute(const MatrixConstView &dataView, Vec1 &steadyState,
|
double compute(const MatrixConstView &dataView, Vec1 &steadyState,
|
||||||
const Mat1 &Q, const Matrix &H, const Vec2 &deepParams,
|
const Mat1 &Q, const Matrix &H, const Vec2 &deepParams,
|
||||||
VectorView &vll, MatrixView &detrendedDataView, size_t start, size_t period,
|
VectorView &vll, MatrixView &detrendedDataView, size_t start, size_t period)
|
||||||
double &penalty, int &info)
|
|
||||||
{
|
{
|
||||||
double lik = INFINITY;
|
double lik = INFINITY;
|
||||||
try
|
|
||||||
{
|
|
||||||
if (period == 0) // initialise all KF matrices
|
if (period == 0) // initialise all KF matrices
|
||||||
initKalmanFilter.initialize(steadyState, deepParams, R, Q, RQRt, T, Pstar, Pinf,
|
initKalmanFilter.initialize(steadyState, deepParams, R, Q, RQRt, T, Pstar, Pinf,
|
||||||
penalty, dataView, detrendedDataView, info);
|
dataView, detrendedDataView);
|
||||||
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, Q, RQRt, T,
|
initKalmanFilter.initialize(steadyState, deepParams, R, Q, RQRt, T,
|
||||||
penalty, dataView, detrendedDataView, info);
|
dataView, detrendedDataView);
|
||||||
|
|
||||||
lik = filter(detrendedDataView, H, vll, start, info);
|
return filter(detrendedDataView, H, vll, start);
|
||||||
}
|
}
|
||||||
catch (const DecisionRules::BlanchardKahnException &bke)
|
|
||||||
{
|
|
||||||
info = 22;
|
|
||||||
return penalty;
|
|
||||||
}
|
|
||||||
|
|
||||||
if (info != 0)
|
|
||||||
return penalty;
|
|
||||||
else
|
|
||||||
return lik;
|
|
||||||
|
|
||||||
};
|
|
||||||
|
|
||||||
private:
|
private:
|
||||||
const std::vector<size_t> zeta_varobs_back_mixed;
|
const std::vector<size_t> zeta_varobs_back_mixed;
|
||||||
|
@ -106,7 +92,7 @@ private:
|
||||||
Vector FUTP; // F upper triangle packed as vector FUTP(i + (j-1)*j/2) = F(i,j) for 1<=i<=j;
|
Vector FUTP; // F upper triangle packed as vector FUTP(i + (j-1)*j/2) = F(i,j) for 1<=i<=j;
|
||||||
|
|
||||||
// Method
|
// Method
|
||||||
double filter(const MatrixView &detrendedDataView, const Matrix &H, VectorView &vll, size_t start, int &info);
|
double filter(const MatrixView &detrendedDataView, const Matrix &H, VectorView &vll, size_t start);
|
||||||
|
|
||||||
};
|
};
|
||||||
|
|
||||||
|
|
|
@ -29,11 +29,11 @@ LogLikelihoodMain::LogLikelihoodMain(const std::string &basename, EstimatedParam
|
||||||
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 double qz_criterium,
|
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> &varobs, double riccati_tol, double lyapunov_tol,
|
const std::vector<size_t> &varobs, double riccati_tol, double lyapunov_tol,
|
||||||
bool noconstant_arg, int &info_arg)
|
bool noconstant_arg)
|
||||||
|
|
||||||
: estSubsamples(estiParDesc.estSubsamples),
|
: estSubsamples(estiParDesc.estSubsamples),
|
||||||
logLikelihoodSubSample(basename, estiParDesc, n_endo, n_exo, zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg, zeta_static_arg, qz_criterium,
|
logLikelihoodSubSample(basename, estiParDesc, n_endo, n_exo, zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg, zeta_static_arg, qz_criterium,
|
||||||
varobs, riccati_tol, lyapunov_tol, noconstant_arg, info_arg),
|
varobs, riccati_tol, lyapunov_tol, noconstant_arg),
|
||||||
vll(estiParDesc.getNumberOfPeriods()), // time dimension size of data
|
vll(estiParDesc.getNumberOfPeriods()), // time dimension size of data
|
||||||
detrendedData(varobs.size(), estiParDesc.getNumberOfPeriods())
|
detrendedData(varobs.size(), estiParDesc.getNumberOfPeriods())
|
||||||
{
|
{
|
||||||
|
|
|
@ -41,7 +41,7 @@ public:
|
||||||
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_arg, const std::vector<size_t> &varobs_arg,
|
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,
|
double riccati_tol_arg, double lyapunov_tol_arg,
|
||||||
bool noconstant_arg, int &info);
|
bool noconstant_arg);
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Compute method Inputs:
|
* Compute method Inputs:
|
||||||
|
@ -54,7 +54,7 @@ public:
|
||||||
|
|
||||||
template <class VEC1, class VEC2>
|
template <class VEC1, class VEC2>
|
||||||
double compute(VEC1 &steadyState, VEC2 &estParams, VectorView &deepParams, const MatrixConstView &data,
|
double compute(VEC1 &steadyState, VEC2 &estParams, VectorView &deepParams, const MatrixConstView &data,
|
||||||
MatrixView &Q, Matrix &H, size_t start, int &info)
|
MatrixView &Q, Matrix &H, size_t start)
|
||||||
{
|
{
|
||||||
double logLikelihood = 0;
|
double logLikelihood = 0;
|
||||||
for (size_t i = 0; i < estSubsamples.size(); ++i)
|
for (size_t i = 0; i < estSubsamples.size(); ++i)
|
||||||
|
@ -68,7 +68,7 @@ public:
|
||||||
|
|
||||||
VectorView vllView(vll, estSubsamples[i].startPeriod, estSubsamples[i].endPeriod-estSubsamples[i].startPeriod+1);
|
VectorView vllView(vll, estSubsamples[i].startPeriod, estSubsamples[i].endPeriod-estSubsamples[i].startPeriod+1);
|
||||||
logLikelihood += logLikelihoodSubSample.compute(vSteadyState, dataView, estParams, deepParams,
|
logLikelihood += logLikelihoodSubSample.compute(vSteadyState, dataView, estParams, deepParams,
|
||||||
Q, H, vllView, detrendedDataView, info, start, i);
|
Q, H, vllView, detrendedDataView, start, i);
|
||||||
}
|
}
|
||||||
return logLikelihood;
|
return logLikelihood;
|
||||||
};
|
};
|
||||||
|
|
|
@ -32,10 +32,10 @@ LogLikelihoodSubSample::~LogLikelihoodSubSample()
|
||||||
LogLikelihoodSubSample::LogLikelihoodSubSample(const std::string &basename, EstimatedParametersDescription &INestiParDesc, size_t n_endo, size_t n_exo,
|
LogLikelihoodSubSample::LogLikelihoodSubSample(const std::string &basename, EstimatedParametersDescription &INestiParDesc, 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 double qz_criterium,
|
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> &varobs, double riccati_tol, double lyapunov_tol, bool noconstant_arg, int &INinfo) :
|
const std::vector<size_t> &varobs, double riccati_tol, double lyapunov_tol, bool noconstant_arg) :
|
||||||
startPenalty(-1e8), estiParDesc(INestiParDesc),
|
estiParDesc(INestiParDesc),
|
||||||
kalmanFilter(basename, n_endo, n_exo, zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg, zeta_static_arg, qz_criterium,
|
kalmanFilter(basename, n_endo, n_exo, zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg, zeta_static_arg, qz_criterium,
|
||||||
varobs, riccati_tol, lyapunov_tol, noconstant_arg, INinfo), eigQ(n_exo), eigH(varobs.size()), info(INinfo)
|
varobs, riccati_tol, lyapunov_tol, noconstant_arg), eigQ(n_exo), eigH(varobs.size())
|
||||||
{
|
{
|
||||||
};
|
};
|
||||||
|
|
||||||
|
|
|
@ -39,35 +39,33 @@ public:
|
||||||
LogLikelihoodSubSample(const std::string &basename, EstimatedParametersDescription &estiParDesc, size_t n_endo, size_t n_exo,
|
LogLikelihoodSubSample(const std::string &basename, 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_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> &zeta_mixed_arg, const std::vector<size_t> &zeta_static_arg, const double qz_criterium,
|
||||||
const std::vector<size_t> &varobs_arg, double riccati_tol_in, double lyapunov_tol, bool noconstant_arg, int &info);
|
const std::vector<size_t> &varobs_arg, double riccati_tol_in, double lyapunov_tol, bool noconstant_arg);
|
||||||
|
|
||||||
template <class VEC1, class VEC2>
|
template <class VEC1, class VEC2>
|
||||||
double compute(VEC1 &steadyState, const MatrixConstView &dataView, VEC2 &estParams, VectorView &deepParams,
|
double compute(VEC1 &steadyState, const MatrixConstView &dataView, VEC2 &estParams, VectorView &deepParams,
|
||||||
MatrixView &Q, Matrix &H, VectorView &vll, MatrixView &detrendedDataView, int &info, size_t start, size_t period)
|
MatrixView &Q, Matrix &H, VectorView &vll, MatrixView &detrendedDataView, size_t start, size_t period)
|
||||||
{
|
{
|
||||||
penalty = startPenalty;
|
|
||||||
logLikelihood = startPenalty;
|
|
||||||
|
|
||||||
updateParams(estParams, deepParams, Q, H, period);
|
updateParams(estParams, deepParams, Q, H, period);
|
||||||
if (info == 0)
|
|
||||||
logLikelihood = kalmanFilter.compute(dataView, steadyState, Q, H, deepParams, vll, detrendedDataView, start, period, penalty, info);
|
|
||||||
// else
|
|
||||||
// logLikelihood+=penalty;
|
|
||||||
|
|
||||||
return logLikelihood;
|
return kalmanFilter.compute(dataView, steadyState, Q, H, deepParams, vll, detrendedDataView, start, period);
|
||||||
|
}
|
||||||
};
|
|
||||||
|
|
||||||
virtual ~LogLikelihoodSubSample();
|
virtual ~LogLikelihoodSubSample();
|
||||||
|
|
||||||
|
class UpdateParamsException
|
||||||
|
{
|
||||||
|
public:
|
||||||
|
double penalty;
|
||||||
|
UpdateParamsException(double penalty_arg) : penalty(penalty_arg)
|
||||||
|
{
|
||||||
|
}
|
||||||
|
};
|
||||||
|
|
||||||
private:
|
private:
|
||||||
double startPenalty, penalty;
|
|
||||||
double logLikelihood;
|
|
||||||
EstimatedParametersDescription &estiParDesc;
|
EstimatedParametersDescription &estiParDesc;
|
||||||
KalmanFilter kalmanFilter;
|
KalmanFilter kalmanFilter;
|
||||||
VDVEigDecomposition eigQ;
|
VDVEigDecomposition eigQ;
|
||||||
VDVEigDecomposition eigH;
|
VDVEigDecomposition eigH;
|
||||||
int &info;
|
|
||||||
|
|
||||||
// methods
|
// methods
|
||||||
template <class VEC>
|
template <class VEC>
|
||||||
|
@ -78,7 +76,6 @@ private:
|
||||||
int test;
|
int test;
|
||||||
bool found;
|
bool found;
|
||||||
std::vector<size_t>::const_iterator it;
|
std::vector<size_t>::const_iterator it;
|
||||||
info = 0;
|
|
||||||
|
|
||||||
for (i = 0; i < estParams.getSize(); ++i)
|
for (i = 0; i < estParams.getSize(); ++i)
|
||||||
{
|
{
|
||||||
|
@ -108,13 +105,9 @@ private:
|
||||||
Q(k2, k1) = Q(k1, k2);
|
Q(k2, k1) = Q(k1, k2);
|
||||||
// [CholQ,testQ] = chol(Q);
|
// [CholQ,testQ] = chol(Q);
|
||||||
test = lapack::choleskyDecomp(Q, "L");
|
test = lapack::choleskyDecomp(Q, "L");
|
||||||
|
assert(test >= 0);
|
||||||
|
|
||||||
if (test > 0)
|
if (test > 0)
|
||||||
{
|
|
||||||
mexPrintf("Caught unhandled exception with cholesky of Q matrix: ");
|
|
||||||
logLikelihood = penalty;
|
|
||||||
info = 1;
|
|
||||||
}
|
|
||||||
else if (test < 0)
|
|
||||||
{
|
{
|
||||||
// The variance-covariance matrix of the structural innovations is not definite positive.
|
// The variance-covariance matrix of the structural innovations is not definite positive.
|
||||||
// We have to compute the eigenvalues of this matrix in order to build the penalty.
|
// We have to compute the eigenvalues of this matrix in order to build the penalty.
|
||||||
|
@ -129,8 +122,7 @@ private:
|
||||||
delta -= evQ(i);
|
delta -= evQ(i);
|
||||||
}
|
}
|
||||||
|
|
||||||
logLikelihood = penalty+delta;
|
throw UpdateParamsException(delta);
|
||||||
info = 43;
|
|
||||||
} // if
|
} // if
|
||||||
break;
|
break;
|
||||||
|
|
||||||
|
@ -144,13 +136,9 @@ private:
|
||||||
|
|
||||||
//[CholH,testH] = chol(H);
|
//[CholH,testH] = chol(H);
|
||||||
test = lapack::choleskyDecomp(H, "L");
|
test = lapack::choleskyDecomp(H, "L");
|
||||||
|
assert(test >= 0);
|
||||||
|
|
||||||
if (test > 0)
|
if (test > 0)
|
||||||
{
|
|
||||||
mexPrintf("Caught unhandled exception with cholesky of Q matrix: ");
|
|
||||||
logLikelihood = penalty;
|
|
||||||
info = 1;
|
|
||||||
}
|
|
||||||
else if (test < 0)
|
|
||||||
{
|
{
|
||||||
// The variance-covariance matrix of the measurement errors is not definite positive.
|
// The variance-covariance matrix of the measurement errors is not definite positive.
|
||||||
// We have to compute the eigenvalues of this matrix in order to build the penalty.
|
// We have to compute the eigenvalues of this matrix in order to build the penalty.
|
||||||
|
@ -165,8 +153,7 @@ private:
|
||||||
if (evH(i) < 0)
|
if (evH(i) < 0)
|
||||||
delta -= evH(i);
|
delta -= evH(i);
|
||||||
}
|
}
|
||||||
logLikelihood = penalty+delta;
|
throw UpdateParamsException(delta);
|
||||||
info = 44;
|
|
||||||
} // end if
|
} // end if
|
||||||
break;
|
break;
|
||||||
|
|
||||||
|
@ -176,8 +163,7 @@ private:
|
||||||
deepParams(k) = estParams(i);
|
deepParams(k) = estParams(i);
|
||||||
break;
|
break;
|
||||||
default:
|
default:
|
||||||
logLikelihood = penalty;
|
assert(false);
|
||||||
info = 1;
|
|
||||||
} // end switch
|
} // end switch
|
||||||
} // end found
|
} // end found
|
||||||
} //end for
|
} //end for
|
||||||
|
|
|
@ -33,10 +33,10 @@ LogPosteriorDensity::LogPosteriorDensity(const std::string &modName, EstimatedPa
|
||||||
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_arg, const std::vector<size_t> &varobs_arg,
|
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,
|
double riccati_tol_arg, double lyapunov_tol_arg,
|
||||||
bool noconstant_arg, int &info_arg) :
|
bool noconstant_arg) :
|
||||||
logPriorDensity(estParamsDesc),
|
logPriorDensity(estParamsDesc),
|
||||||
logLikelihoodMain(modName, estParamsDesc, n_endo, n_exo, zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg,
|
logLikelihoodMain(modName, estParamsDesc, 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, noconstant_arg, info_arg)
|
zeta_static_arg, qz_criterium_arg, varobs_arg, riccati_tol_arg, lyapunov_tol_arg, noconstant_arg)
|
||||||
{
|
{
|
||||||
|
|
||||||
}
|
}
|
||||||
|
|
|
@ -48,13 +48,13 @@ public:
|
||||||
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_arg, const std::vector<size_t> &varobs_arg,
|
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,
|
double riccati_tol_arg, double lyapunov_tol_arg,
|
||||||
bool noconstant_arg, int &info_arg);
|
bool noconstant_arg);
|
||||||
|
|
||||||
template <class VEC1, class VEC2>
|
template <class VEC1, class VEC2>
|
||||||
double
|
double
|
||||||
compute(VEC1 &steadyState, VEC2 &estParams, VectorView &deepParams, const MatrixConstView &data, MatrixView &Q, Matrix &H, size_t presampleStart, int &info)
|
compute(VEC1 &steadyState, VEC2 &estParams, VectorView &deepParams, const MatrixConstView &data, MatrixView &Q, Matrix &H, size_t presampleStart)
|
||||||
{
|
{
|
||||||
return -logLikelihoodMain.compute(steadyState, estParams, deepParams, data, Q, H, presampleStart, info)
|
return -logLikelihoodMain.compute(steadyState, estParams, deepParams, data, Q, H, presampleStart)
|
||||||
-logPriorDensity.compute(estParams);
|
-logPriorDensity.compute(estParams);
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
|
@ -1,5 +1,5 @@
|
||||||
/*
|
/*
|
||||||
* Copyright (C) 2010-2012 Dynare Team
|
* Copyright (C) 2010-2013 Dynare Team
|
||||||
*
|
*
|
||||||
* This file is part of Dynare.
|
* This file is part of Dynare.
|
||||||
*
|
*
|
||||||
|
@ -47,7 +47,7 @@ public:
|
||||||
template<class VEC1>
|
template<class VEC1>
|
||||||
double compute(VectorView &mhLogPostDens, MatrixView &mhParams, VEC1 &steadyState,
|
double compute(VectorView &mhLogPostDens, MatrixView &mhParams, VEC1 &steadyState,
|
||||||
Vector &estParams, VectorView &deepParams, const MatrixConstView &data, MatrixView &Q, Matrix &H,
|
Vector &estParams, VectorView &deepParams, const MatrixConstView &data, MatrixView &Q, Matrix &H,
|
||||||
const size_t presampleStart, int &info, const size_t startDraw, size_t nMHruns,
|
const size_t presampleStart, const size_t startDraw, size_t nMHruns,
|
||||||
LogPosteriorDensity &lpd, Proposal &pDD, EstimatedParametersDescription &epd)
|
LogPosteriorDensity &lpd, Proposal &pDD, EstimatedParametersDescription &epd)
|
||||||
{
|
{
|
||||||
//streambuf *likbuf, *drawbuf *backup;
|
//streambuf *likbuf, *drawbuf *backup;
|
||||||
|
@ -60,7 +60,7 @@ public:
|
||||||
size_t count, accepted = 0;
|
size_t count, accepted = 0;
|
||||||
parDraw = estParams;
|
parDraw = estParams;
|
||||||
|
|
||||||
logpost = -lpd.compute(steadyState, estParams, deepParams, data, Q, H, presampleStart, info);
|
logpost = -lpd.compute(steadyState, estParams, deepParams, data, Q, H, presampleStart);
|
||||||
|
|
||||||
for (size_t run = startDraw - 1; run < nMHruns; ++run)
|
for (size_t run = startDraw - 1; run < nMHruns; ++run)
|
||||||
{
|
{
|
||||||
|
@ -79,7 +79,7 @@ public:
|
||||||
{
|
{
|
||||||
try
|
try
|
||||||
{
|
{
|
||||||
newLogpost = -lpd.compute(steadyState, newParDraw, deepParams, data, Q, H, presampleStart, info);
|
newLogpost = -lpd.compute(steadyState, newParDraw, deepParams, data, Q, H, presampleStart);
|
||||||
}
|
}
|
||||||
catch (const std::exception &e)
|
catch (const std::exception &e)
|
||||||
{
|
{
|
||||||
|
|
|
@ -114,7 +114,7 @@ fillEstParamsInfo(const mxArray *estim_params_info, EstimatedParameter::pType ty
|
||||||
int
|
int
|
||||||
sampleMHMC(LogPosteriorDensity &lpd, RandomWalkMetropolisHastings &rwmh,
|
sampleMHMC(LogPosteriorDensity &lpd, RandomWalkMetropolisHastings &rwmh,
|
||||||
VectorView &steadyState, VectorConstView &estParams, VectorView &deepParams, const MatrixConstView &data,
|
VectorView &steadyState, VectorConstView &estParams, VectorView &deepParams, const MatrixConstView &data,
|
||||||
MatrixView &Q, Matrix &H, size_t presampleStart, int &info, const VectorConstView &nruns,
|
MatrixView &Q, Matrix &H, size_t presampleStart, const VectorConstView &nruns,
|
||||||
size_t fblock, size_t nBlocks, Proposal pdd, EstimatedParametersDescription &epd,
|
size_t fblock, size_t nBlocks, Proposal pdd, EstimatedParametersDescription &epd,
|
||||||
const std::string &resultsFileStem, size_t console_mode, size_t load_mh_file)
|
const std::string &resultsFileStem, size_t console_mode, size_t load_mh_file)
|
||||||
{
|
{
|
||||||
|
@ -397,7 +397,7 @@ sampleMHMC(LogPosteriorDensity &lpd, RandomWalkMetropolisHastings &rwmh,
|
||||||
try
|
try
|
||||||
{
|
{
|
||||||
jsux = rwmh.compute(mhLogPostDens, mhParamDraws, steadyState, startParams, deepParams, data, Q, H,
|
jsux = rwmh.compute(mhLogPostDens, mhParamDraws, steadyState, startParams, deepParams, data, Q, H,
|
||||||
presampleStart, info, irun, currInitSizeArray, lpd, pdd, epd);
|
presampleStart, irun, currInitSizeArray, lpd, pdd, epd);
|
||||||
irun = currInitSizeArray;
|
irun = currInitSizeArray;
|
||||||
sux += jsux*currInitSizeArray;
|
sux += jsux*currInitSizeArray;
|
||||||
j += currInitSizeArray; //j=j+1;
|
j += currInitSizeArray; //j=j+1;
|
||||||
|
@ -715,7 +715,7 @@ logMCMCposterior(VectorConstView &estParams, const MatrixConstView &data,
|
||||||
// Allocate LogPosteriorDensity object
|
// Allocate LogPosteriorDensity object
|
||||||
int info;
|
int info;
|
||||||
LogPosteriorDensity lpd(basename, epd, n_endo, n_exo, zeta_fwrd, zeta_back, zeta_mixed, zeta_static,
|
LogPosteriorDensity lpd(basename, epd, n_endo, n_exo, zeta_fwrd, zeta_back, zeta_mixed, zeta_static,
|
||||||
qz_criterium, varobs, riccati_tol, lyapunov_tol, noconstant, info);
|
qz_criterium, varobs, riccati_tol, lyapunov_tol, noconstant);
|
||||||
|
|
||||||
// Construct MHMCMC Sampler
|
// Construct MHMCMC Sampler
|
||||||
RandomWalkMetropolisHastings rwmh(estParams.getSize());
|
RandomWalkMetropolisHastings rwmh(estParams.getSize());
|
||||||
|
@ -727,7 +727,7 @@ logMCMCposterior(VectorConstView &estParams, const MatrixConstView &data,
|
||||||
Proposal pdd(vJscale, D);
|
Proposal pdd(vJscale, D);
|
||||||
|
|
||||||
//sample MHMCMC draws and get get last line run in the last MH block sub-array
|
//sample MHMCMC draws and get get last line run in the last MH block sub-array
|
||||||
int lastMHblockArrayLine = sampleMHMC(lpd, rwmh, steadyState, estParams, deepParams, data, Q, H, presample, info,
|
int lastMHblockArrayLine = sampleMHMC(lpd, rwmh, steadyState, estParams, deepParams, data, Q, H, presample,
|
||||||
nMHruns, fblock, nBlocks, pdd, epd, resultsFileStem, console_mode, load_mh_file);
|
nMHruns, fblock, nBlocks, pdd, epd, resultsFileStem, console_mode, load_mh_file);
|
||||||
|
|
||||||
// Cleanups
|
// Cleanups
|
||||||
|
|
|
@ -104,7 +104,7 @@ double
|
||||||
logposterior(VEC1 &estParams, const MatrixConstView &data,
|
logposterior(VEC1 &estParams, const MatrixConstView &data,
|
||||||
const mxArray *options_, const mxArray *M_, const mxArray *estim_params_,
|
const mxArray *options_, const mxArray *M_, const mxArray *estim_params_,
|
||||||
const mxArray *bayestopt_, const mxArray *oo_, VEC2 &steadyState, double *trend_coeff,
|
const mxArray *bayestopt_, const mxArray *oo_, VEC2 &steadyState, double *trend_coeff,
|
||||||
int &info, VectorView &deepParams, Matrix &H, MatrixView &Q)
|
VectorView &deepParams, Matrix &H, MatrixView &Q)
|
||||||
{
|
{
|
||||||
double loglinear = *mxGetPr(mxGetField(options_, 0, "loglinear"));
|
double loglinear = *mxGetPr(mxGetField(options_, 0, "loglinear"));
|
||||||
if (loglinear == 1)
|
if (loglinear == 1)
|
||||||
|
@ -183,12 +183,12 @@ logposterior(VEC1 &estParams, const MatrixConstView &data,
|
||||||
|
|
||||||
// Allocate LogPosteriorDensity object
|
// Allocate LogPosteriorDensity object
|
||||||
LogPosteriorDensity lpd(basename, epd, n_endo, n_exo, zeta_fwrd, zeta_back, zeta_mixed, zeta_static,
|
LogPosteriorDensity lpd(basename, epd, n_endo, n_exo, zeta_fwrd, zeta_back, zeta_mixed, zeta_static,
|
||||||
qz_criterium, varobs, riccati_tol, lyapunov_tol, noconstant, info);
|
qz_criterium, varobs, riccati_tol, lyapunov_tol, noconstant);
|
||||||
|
|
||||||
// Construct arguments of compute() method
|
// Construct arguments of compute() method
|
||||||
|
|
||||||
// Compute the posterior
|
// Compute the posterior
|
||||||
double logPD = lpd.compute(steadyState, estParams, deepParams, data, Q, H, presample, info);
|
double logPD = lpd.compute(steadyState, estParams, deepParams, data, Q, H, presample);
|
||||||
|
|
||||||
// Cleanups
|
// Cleanups
|
||||||
for (std::vector<EstimatedParameter>::iterator it = estParamsInfo.begin();
|
for (std::vector<EstimatedParameter>::iterator it = estParamsInfo.begin();
|
||||||
|
@ -268,11 +268,10 @@ mexFunction(int nlhs, mxArray *plhs[],
|
||||||
// Compute and return the value
|
// Compute and return the value
|
||||||
try
|
try
|
||||||
{
|
{
|
||||||
int info;
|
|
||||||
*lik = logposterior(estParams, data, options_, M_, estim_params_, bayestopt_, oo_,
|
*lik = logposterior(estParams, data, options_, M_, estim_params_, bayestopt_, oo_,
|
||||||
steadyState, trend_coeff, info, deepParams, H, Q);
|
steadyState, trend_coeff, deepParams, H, Q);
|
||||||
*info_mx = info;
|
*info_mx = 0;
|
||||||
*exit_flag = info;
|
*exit_flag = 0;
|
||||||
}
|
}
|
||||||
catch (LogposteriorMexErrMsgTxtException e)
|
catch (LogposteriorMexErrMsgTxtException e)
|
||||||
{
|
{
|
||||||
|
|
Loading…
Reference in New Issue