updating calling sequence for estimation dll and introducing templates
when necessarytime-shift
parent
eb567f5202
commit
da0beac147
|
@ -63,6 +63,5 @@ nodist_logMHMCMCposterior_SOURCES = \
|
||||||
$(COMMON_SRCS) \
|
$(COMMON_SRCS) \
|
||||||
$(TOPDIR)/Proposal.cc \
|
$(TOPDIR)/Proposal.cc \
|
||||||
$(TOPDIR)/Proposal.hh \
|
$(TOPDIR)/Proposal.hh \
|
||||||
$(TOPDIR)/RandomWalkMetropolisHastings.cc \
|
|
||||||
$(TOPDIR)/RandomWalkMetropolisHastings.hh \
|
$(TOPDIR)/RandomWalkMetropolisHastings.hh \
|
||||||
$(TOPDIR)/logMHMCMCposterior.cc
|
$(TOPDIR)/logMHMCMCposterior.cc
|
||||||
|
|
|
@ -55,31 +55,6 @@ InitializeKalmanFilter::InitializeKalmanFilter(const std::string &dynamicDllFile
|
||||||
zeta_back_mixed[i]) - zeta_varobs_back_mixed.begin());
|
zeta_back_mixed[i]) - zeta_varobs_back_mixed.begin());
|
||||||
|
|
||||||
}
|
}
|
||||||
// initialise parameter dependent KF matrices only but not Ps
|
|
||||||
void
|
|
||||||
InitializeKalmanFilter::initialize(VectorView &steadyState, const Vector &deepParams, Matrix &R,
|
|
||||||
const Matrix &Q, Matrix &RQRt, Matrix &T,
|
|
||||||
double &penalty, const MatrixConstView &dataView,
|
|
||||||
MatrixView &detrendedDataView, int &info)
|
|
||||||
{
|
|
||||||
modelSolution.compute(steadyState, deepParams, g_x, g_u);
|
|
||||||
detrendData.detrend(steadyState, dataView, detrendedDataView);
|
|
||||||
|
|
||||||
setT(T, info);
|
|
||||||
setRQR(R, Q, RQRt, info);
|
|
||||||
}
|
|
||||||
|
|
||||||
// initialise all KF matrices
|
|
||||||
void
|
|
||||||
InitializeKalmanFilter::initialize(VectorView &steadyState, const Vector &deepParams, Matrix &R,
|
|
||||||
const Matrix &Q, Matrix &RQRt, Matrix &T, Matrix &Pstar, Matrix &Pinf,
|
|
||||||
double &penalty, const MatrixConstView &dataView,
|
|
||||||
MatrixView &detrendedDataView, int &info)
|
|
||||||
{
|
|
||||||
initialize(steadyState, deepParams, R, Q, RQRt, T, penalty, dataView, detrendedDataView, info);
|
|
||||||
setPstar(Pstar, Pinf, T, RQRt, info);
|
|
||||||
}
|
|
||||||
|
|
||||||
void
|
void
|
||||||
InitializeKalmanFilter::setT(Matrix &T, int &info)
|
InitializeKalmanFilter::setT(Matrix &T, int &info)
|
||||||
{
|
{
|
||||||
|
@ -89,7 +64,7 @@ InitializeKalmanFilter::setT(Matrix &T, int &info)
|
||||||
}
|
}
|
||||||
|
|
||||||
void
|
void
|
||||||
InitializeKalmanFilter::setRQR(Matrix &R, const Matrix &Q, Matrix &RQRt, int &info)
|
InitializeKalmanFilter::setRQR(Matrix &R, const MatrixView &Q, Matrix &RQRt, int &info)
|
||||||
{
|
{
|
||||||
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);
|
||||||
|
|
||||||
|
|
|
@ -48,12 +48,30 @@ public:
|
||||||
const std::vector<size_t> &zeta_varobs_back_mixed_arg,
|
const std::vector<size_t> &zeta_varobs_back_mixed_arg,
|
||||||
double qz_criterium_arg, double lyapunov_tol_arg, int &info);
|
double qz_criterium_arg, double lyapunov_tol_arg, int &info);
|
||||||
virtual ~InitializeKalmanFilter();
|
virtual ~InitializeKalmanFilter();
|
||||||
// initialise all KF matrices
|
|
||||||
void initialize(VectorView &steadyState, const Vector &deepParams, Matrix &R, const Matrix &Q, Matrix &RQRt,
|
|
||||||
Matrix &T, Matrix &Pstar, Matrix &Pinf, double &penalty, const MatrixConstView &dataView, MatrixView &detrendedDataView, int &info);
|
|
||||||
// initialise parameter dependent KF matrices only but not Ps
|
// initialise parameter dependent KF matrices only but not Ps
|
||||||
void initialize(VectorView &steadyState, const Vector &deepParams, Matrix &R, const Matrix &Q, Matrix &RQRt,
|
template <class VEC>
|
||||||
Matrix &T, double &penalty, const MatrixConstView &dataView, MatrixView &detrendedDataView, int &info);
|
void initialize(VEC &steadyState, const VectorView &deepParams, Matrix &R,
|
||||||
|
const MatrixView &Q, Matrix &RQRt, Matrix &T,
|
||||||
|
double &penalty, const MatrixConstView &dataView,
|
||||||
|
MatrixView &detrendedDataView, int &info)
|
||||||
|
{
|
||||||
|
modelSolution.compute(steadyState, deepParams, g_x, g_u);
|
||||||
|
detrendData.detrend(steadyState, dataView, detrendedDataView);
|
||||||
|
|
||||||
|
setT(T, info);
|
||||||
|
setRQR(R, Q, RQRt, info);
|
||||||
|
}
|
||||||
|
|
||||||
|
// initialise all KF matrices
|
||||||
|
template <class VEC>
|
||||||
|
void initialize(VEC &steadyState, const VectorView &deepParams, Matrix &R,
|
||||||
|
const MatrixView &Q, Matrix &RQRt, Matrix &T, Matrix &Pstar, Matrix &Pinf,
|
||||||
|
double &penalty, const MatrixConstView &dataView,
|
||||||
|
MatrixView &detrendedDataView, int &info)
|
||||||
|
{
|
||||||
|
initialize(steadyState, deepParams, R, Q, RQRt, T, penalty, dataView, detrendedDataView, info);
|
||||||
|
setPstar(Pstar, Pinf, T, RQRt, info);
|
||||||
|
}
|
||||||
|
|
||||||
private:
|
private:
|
||||||
const double lyapunov_tol;
|
const double lyapunov_tol;
|
||||||
|
@ -68,7 +86,7 @@ private:
|
||||||
Matrix g_u;
|
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 MatrixView &Q, Matrix &RQRt, int &info);
|
||||||
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, int &info);
|
||||||
|
|
||||||
};
|
};
|
||||||
|
|
|
@ -73,36 +73,6 @@ KalmanFilter::compute_zeta_varobs_back_mixed(const std::vector<size_t> &zeta_bac
|
||||||
return zeta_varobs_back_mixed;
|
return zeta_varobs_back_mixed;
|
||||||
}
|
}
|
||||||
|
|
||||||
double
|
|
||||||
KalmanFilter::compute(const MatrixConstView &dataView, VectorView &steadyState,
|
|
||||||
const Matrix &Q, const Matrix &H, const Vector &deepParams,
|
|
||||||
VectorView &vll, MatrixView &detrendedDataView,
|
|
||||||
size_t start, size_t period, double &penalty, int &info)
|
|
||||||
{
|
|
||||||
double lik = INFINITY;
|
|
||||||
try
|
|
||||||
{
|
|
||||||
if (period == 0) // initialise all KF matrices
|
|
||||||
initKalmanFilter.initialize(steadyState, deepParams, R, Q, RQRt, T, Pstar, Pinf,
|
|
||||||
penalty, dataView, detrendedDataView, info);
|
|
||||||
else // initialise parameter dependent KF matrices only but not Ps
|
|
||||||
initKalmanFilter.initialize(steadyState, deepParams, R, Q, RQRt, T,
|
|
||||||
penalty, dataView, detrendedDataView, info);
|
|
||||||
|
|
||||||
lik = filter(detrendedDataView, H, vll, start, info);
|
|
||||||
}
|
|
||||||
catch (const DecisionRules::BlanchardKahnException &bke)
|
|
||||||
{
|
|
||||||
info = 22;
|
|
||||||
return penalty;
|
|
||||||
}
|
|
||||||
|
|
||||||
if (info != 0)
|
|
||||||
return penalty;
|
|
||||||
else
|
|
||||||
return lik;
|
|
||||||
|
|
||||||
};
|
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Multi-variate standard Kalman Filter
|
* Multi-variate standard Kalman Filter
|
||||||
|
@ -110,7 +80,7 @@ KalmanFilter::compute(const MatrixConstView &dataView, VectorView &steadyState,
|
||||||
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, int &info)
|
||||||
{
|
{
|
||||||
double loglik = 0.0, ll, logFdet, 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);
|
||||||
|
|
|
@ -54,10 +54,36 @@ public:
|
||||||
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, int &info);
|
double riccati_tol_arg, double lyapunov_tol_arg, int &info);
|
||||||
|
|
||||||
double compute(const MatrixConstView &dataView, VectorView &steadyState,
|
template <class VEC>
|
||||||
const Matrix &Q, const Matrix &H, const Vector &deepParams,
|
double compute(const MatrixConstView &dataView, VEC &steadyState,
|
||||||
|
const MatrixView &Q, const Matrix &H, const VectorView &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 &penalty, int &info)
|
||||||
|
{
|
||||||
|
double lik = INFINITY;
|
||||||
|
try
|
||||||
|
{
|
||||||
|
if (period == 0) // initialise all KF matrices
|
||||||
|
initKalmanFilter.initialize(steadyState, deepParams, R, Q, RQRt, T, Pstar, Pinf,
|
||||||
|
penalty, dataView, detrendedDataView, info);
|
||||||
|
else // initialise parameter dependent KF matrices only but not Ps
|
||||||
|
initKalmanFilter.initialize(steadyState, deepParams, R, Q, RQRt, T,
|
||||||
|
penalty, dataView, detrendedDataView, info);
|
||||||
|
|
||||||
|
lik = filter(detrendedDataView, H, vll, start, info);
|
||||||
|
}
|
||||||
|
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;
|
||||||
|
|
|
@ -44,23 +44,4 @@ LogLikelihoodMain::~LogLikelihoodMain()
|
||||||
|
|
||||||
}
|
}
|
||||||
|
|
||||||
double
|
|
||||||
LogLikelihoodMain::compute(Matrix &steadyState, const Vector &estParams, Vector &deepParams, const MatrixConstView &data, Matrix &Q, Matrix &H, size_t start, int &info)
|
|
||||||
{
|
|
||||||
double logLikelihood = 0;
|
|
||||||
for (size_t i = 0; i < estSubsamples.size(); ++i)
|
|
||||||
{
|
|
||||||
VectorView vSteadyState = mat::get_col(steadyState, i);
|
|
||||||
|
|
||||||
MatrixConstView dataView(data, 0, estSubsamples[i].startPeriod,
|
|
||||||
data.getRows(), estSubsamples[i].endPeriod-estSubsamples[i].startPeriod+1);
|
|
||||||
MatrixView detrendedDataView(detrendedData, 0, estSubsamples[i].startPeriod,
|
|
||||||
data.getRows(), 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,
|
|
||||||
Q, H, vllView, detrendedDataView, info, start, i);
|
|
||||||
}
|
|
||||||
return logLikelihood;
|
|
||||||
};
|
|
||||||
|
|
||||||
|
|
|
@ -50,14 +50,29 @@ public:
|
||||||
* Q and H KF matrices of shock and measurement error varinaces and covariances
|
* Q and H KF matrices of shock and measurement error varinaces and covariances
|
||||||
* KF logLikelihood calculation start period.
|
* KF logLikelihood calculation start period.
|
||||||
*/
|
*/
|
||||||
double compute(Matrix &steadyState, const Vector &estParams, Vector &deepParams, const MatrixConstView &data,
|
|
||||||
Matrix &Q, Matrix &H, size_t presampleStart, int &info); // for calls from estimation and to set Steady State
|
|
||||||
|
|
||||||
Vector &
|
template <class VEC1, class VEC2>
|
||||||
getVll()
|
double compute(VEC1 &steadyState, VEC2 &estParams, VectorView &deepParams, const MatrixConstView &data,
|
||||||
|
MatrixView &Q, Matrix &H, size_t start, int &info)
|
||||||
{
|
{
|
||||||
return vll;
|
double logLikelihood = 0;
|
||||||
|
for (size_t i = 0; i < estSubsamples.size(); ++i)
|
||||||
|
{
|
||||||
|
VectorView vSteadyState (steadyState,0,steadyState.getSize());
|
||||||
|
|
||||||
|
MatrixConstView dataView(data, 0, estSubsamples[i].startPeriod,
|
||||||
|
data.getRows(), estSubsamples[i].endPeriod-estSubsamples[i].startPeriod+1);
|
||||||
|
MatrixView detrendedDataView(detrendedData, 0, estSubsamples[i].startPeriod,
|
||||||
|
data.getRows(), 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,
|
||||||
|
Q, H, vllView, detrendedDataView, info, start, i);
|
||||||
|
}
|
||||||
|
return logLikelihood;
|
||||||
};
|
};
|
||||||
|
|
||||||
|
Vector &getVll() { return vll; };
|
||||||
};
|
};
|
||||||
|
|
||||||
#endif // !defined(E126AEF5_AC28_400a_821A_3BCFD1BC4C22__INCLUDED_)
|
#endif // !defined(E126AEF5_AC28_400a_821A_3BCFD1BC4C22__INCLUDED_)
|
||||||
|
|
|
@ -25,8 +25,6 @@
|
||||||
|
|
||||||
//#include "LogLikelihoodSubSample.hh"
|
//#include "LogLikelihoodSubSample.hh"
|
||||||
#include "LogLikelihoodMain.hh" // use ...Main.hh for testing only
|
#include "LogLikelihoodMain.hh" // use ...Main.hh for testing only
|
||||||
#include <algorithm>
|
|
||||||
#include "LapackBindings.hh"
|
|
||||||
|
|
||||||
LogLikelihoodSubSample::~LogLikelihoodSubSample()
|
LogLikelihoodSubSample::~LogLikelihoodSubSample()
|
||||||
{
|
{
|
||||||
|
@ -42,133 +40,3 @@ LogLikelihoodSubSample::LogLikelihoodSubSample(const std::string &dynamicDllFile
|
||||||
{
|
{
|
||||||
};
|
};
|
||||||
|
|
||||||
double
|
|
||||||
LogLikelihoodSubSample::compute(VectorView &steadyState, const MatrixConstView &dataView, const Vector &estParams, Vector &deepParams,
|
|
||||||
Matrix &Q, Matrix &H, VectorView &vll, MatrixView &detrendedDataView, int &info, size_t start, size_t period)
|
|
||||||
{
|
|
||||||
penalty = startPenalty;
|
|
||||||
logLikelihood = startPenalty;
|
|
||||||
|
|
||||||
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;
|
|
||||||
|
|
||||||
};
|
|
||||||
|
|
||||||
void
|
|
||||||
LogLikelihoodSubSample::updateParams(const Vector &estParams, Vector &deepParams,
|
|
||||||
Matrix &Q, Matrix &H, size_t period)
|
|
||||||
{
|
|
||||||
size_t i, k, k1, k2;
|
|
||||||
int test;
|
|
||||||
bool found;
|
|
||||||
std::vector<size_t>::const_iterator it;
|
|
||||||
info = 0;
|
|
||||||
|
|
||||||
for (i = 0; i < estParams.getSize(); ++i)
|
|
||||||
{
|
|
||||||
found = false;
|
|
||||||
it = find(estiParDesc.estParams[i].subSampleIDs.begin(),
|
|
||||||
estiParDesc.estParams[i].subSampleIDs.end(), period);
|
|
||||||
if (it != estiParDesc.estParams[i].subSampleIDs.end())
|
|
||||||
found = true;
|
|
||||||
if (found)
|
|
||||||
{
|
|
||||||
switch (estiParDesc.estParams[i].ptype)
|
|
||||||
{
|
|
||||||
case EstimatedParameter::shock_SD:
|
|
||||||
k = estiParDesc.estParams[i].ID1;
|
|
||||||
Q(k, k) = estParams(i)*estParams(i);
|
|
||||||
break;
|
|
||||||
|
|
||||||
case EstimatedParameter::measureErr_SD:
|
|
||||||
k = estiParDesc.estParams[i].ID1;
|
|
||||||
H(k, k) = estParams(i)*estParams(i);
|
|
||||||
break;
|
|
||||||
|
|
||||||
case EstimatedParameter::shock_Corr:
|
|
||||||
k1 = estiParDesc.estParams[i].ID1;
|
|
||||||
k2 = estiParDesc.estParams[i].ID2;
|
|
||||||
Q(k1, k2) = estParams(i)*sqrt(Q(k1, k1)*Q(k2, k2));
|
|
||||||
Q(k2, k1) = Q(k1, k2);
|
|
||||||
// [CholQ,testQ] = chol(Q);
|
|
||||||
test = lapack::choleskyDecomp(Q, "L");
|
|
||||||
if (test > 0)
|
|
||||||
{
|
|
||||||
mexPrintf("Caugth 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.
|
|
||||||
// We have to compute the eigenvalues of this matrix in order to build the penalty.
|
|
||||||
double delta = 0;
|
|
||||||
eigQ.calculate(Q); // get eigenvalues
|
|
||||||
//k = find(a < 0);
|
|
||||||
if (eigQ.hasConverged())
|
|
||||||
{
|
|
||||||
const Vector &evQ = eigQ.getD();
|
|
||||||
for (i = 0; i < evQ.getSize(); ++i)
|
|
||||||
if (evQ(i) < 0)
|
|
||||||
delta -= evQ(i);
|
|
||||||
}
|
|
||||||
|
|
||||||
logLikelihood = penalty+delta;
|
|
||||||
info = 43;
|
|
||||||
} // if
|
|
||||||
break;
|
|
||||||
|
|
||||||
case EstimatedParameter::measureErr_Corr:
|
|
||||||
k1 = estiParDesc.estParams[i].ID1;
|
|
||||||
k2 = estiParDesc.estParams[i].ID2;
|
|
||||||
// H(k1,k2) = xparam1(i)*sqrt(H(k1,k1)*H(k2,k2));
|
|
||||||
// H(k2,k1) = H(k1,k2);
|
|
||||||
H(k1, k2) = estParams(i)*sqrt(H(k1, k1)*H(k2, k2));
|
|
||||||
H(k2, k1) = H(k1, k2);
|
|
||||||
|
|
||||||
//[CholH,testH] = chol(H);
|
|
||||||
test = lapack::choleskyDecomp(H, "L");
|
|
||||||
if (test > 0)
|
|
||||||
{
|
|
||||||
mexPrintf("Caugth 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.
|
|
||||||
// We have to compute the eigenvalues of this matrix in order to build the penalty.
|
|
||||||
//a = diag(eig(H));
|
|
||||||
double delta = 0;
|
|
||||||
eigH.calculate(H); // get eigenvalues
|
|
||||||
//k = find(a < 0);
|
|
||||||
if (eigH.hasConverged())
|
|
||||||
{
|
|
||||||
const Vector &evH = eigH.getD();
|
|
||||||
for (i = 0; i < evH.getSize(); ++i)
|
|
||||||
if (evH(i) < 0)
|
|
||||||
delta -= evH(i);
|
|
||||||
}
|
|
||||||
logLikelihood = penalty+delta;
|
|
||||||
info = 44;
|
|
||||||
} // end if
|
|
||||||
break;
|
|
||||||
|
|
||||||
//if estim_params_.np > 0 // i.e. num of deep parameters >0
|
|
||||||
case EstimatedParameter::deepPar:
|
|
||||||
k = estiParDesc.estParams[i].ID1;
|
|
||||||
deepParams(k) = estParams(i);
|
|
||||||
break;
|
|
||||||
default:
|
|
||||||
logLikelihood = penalty;
|
|
||||||
info = 1;
|
|
||||||
} // end switch
|
|
||||||
} // end found
|
|
||||||
} //end for
|
|
||||||
};
|
|
||||||
|
|
||||||
|
|
|
@ -26,9 +26,11 @@
|
||||||
#if !defined(DF8B7AF5_8169_4587_9037_2CD2C82E2DDF__INCLUDED_)
|
#if !defined(DF8B7AF5_8169_4587_9037_2CD2C82E2DDF__INCLUDED_)
|
||||||
#define DF8B7AF5_8169_4587_9037_2CD2C82E2DDF__INCLUDED_
|
#define DF8B7AF5_8169_4587_9037_2CD2C82E2DDF__INCLUDED_
|
||||||
|
|
||||||
|
#include <algorithm>
|
||||||
#include "EstimatedParametersDescription.hh"
|
#include "EstimatedParametersDescription.hh"
|
||||||
#include "KalmanFilter.hh"
|
#include "KalmanFilter.hh"
|
||||||
#include "VDVEigDecomposition.hh"
|
#include "VDVEigDecomposition.hh"
|
||||||
|
#include "LapackBindings.hh"
|
||||||
|
|
||||||
class LogLikelihoodSubSample
|
class LogLikelihoodSubSample
|
||||||
{
|
{
|
||||||
|
@ -39,8 +41,23 @@ public:
|
||||||
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, int &info);
|
const std::vector<size_t> &varobs_arg, double riccati_tol_in, double lyapunov_tol, int &info);
|
||||||
|
|
||||||
double compute(VectorView &steadyState, const MatrixConstView &dataView, const Vector &estParams, Vector &deepParams,
|
template <class VEC1, class VEC2>
|
||||||
Matrix &Q, Matrix &H, VectorView &vll, MatrixView &detrendedDataView, int &info, size_t start, size_t period);
|
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)
|
||||||
|
{
|
||||||
|
penalty = startPenalty;
|
||||||
|
logLikelihood = startPenalty;
|
||||||
|
|
||||||
|
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;
|
||||||
|
|
||||||
|
};
|
||||||
|
|
||||||
virtual ~LogLikelihoodSubSample();
|
virtual ~LogLikelihoodSubSample();
|
||||||
|
|
||||||
private:
|
private:
|
||||||
|
@ -53,8 +70,119 @@ private:
|
||||||
int &info;
|
int &info;
|
||||||
|
|
||||||
// methods
|
// methods
|
||||||
void updateParams(const Vector &estParams, Vector &deepParams,
|
template <class VEC>
|
||||||
Matrix &Q, Matrix &H, size_t period);
|
void updateParams(VEC &estParams, VectorView &deepParams,
|
||||||
|
MatrixView &Q, Matrix &H, size_t period)
|
||||||
|
{
|
||||||
|
size_t i, k, k1, k2;
|
||||||
|
int test;
|
||||||
|
bool found;
|
||||||
|
std::vector<size_t>::const_iterator it;
|
||||||
|
info = 0;
|
||||||
|
|
||||||
|
for (i = 0; i < estParams.getSize(); ++i)
|
||||||
|
{
|
||||||
|
found = false;
|
||||||
|
it = find(estiParDesc.estParams[i].subSampleIDs.begin(),
|
||||||
|
estiParDesc.estParams[i].subSampleIDs.end(), period);
|
||||||
|
if (it != estiParDesc.estParams[i].subSampleIDs.end())
|
||||||
|
found = true;
|
||||||
|
if (found)
|
||||||
|
{
|
||||||
|
switch (estiParDesc.estParams[i].ptype)
|
||||||
|
{
|
||||||
|
case EstimatedParameter::shock_SD:
|
||||||
|
k = estiParDesc.estParams[i].ID1;
|
||||||
|
Q(k, k) = estParams(i)*estParams(i);
|
||||||
|
break;
|
||||||
|
|
||||||
|
case EstimatedParameter::measureErr_SD:
|
||||||
|
k = estiParDesc.estParams[i].ID1;
|
||||||
|
H(k, k) = estParams(i)*estParams(i);
|
||||||
|
break;
|
||||||
|
|
||||||
|
case EstimatedParameter::shock_Corr:
|
||||||
|
k1 = estiParDesc.estParams[i].ID1;
|
||||||
|
k2 = estiParDesc.estParams[i].ID2;
|
||||||
|
Q(k1, k2) = estParams(i)*sqrt(Q(k1, k1)*Q(k2, k2));
|
||||||
|
Q(k2, k1) = Q(k1, k2);
|
||||||
|
// [CholQ,testQ] = chol(Q);
|
||||||
|
test = lapack::choleskyDecomp(Q, "L");
|
||||||
|
if (test > 0)
|
||||||
|
{
|
||||||
|
mexPrintf("Caugth 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.
|
||||||
|
// We have to compute the eigenvalues of this matrix in order to build the penalty.
|
||||||
|
double delta = 0;
|
||||||
|
eigQ.calculate(Q); // get eigenvalues
|
||||||
|
//k = find(a < 0);
|
||||||
|
if (eigQ.hasConverged())
|
||||||
|
{
|
||||||
|
const Vector &evQ = eigQ.getD();
|
||||||
|
for (i = 0; i < evQ.getSize(); ++i)
|
||||||
|
if (evQ(i) < 0)
|
||||||
|
delta -= evQ(i);
|
||||||
|
}
|
||||||
|
|
||||||
|
logLikelihood = penalty+delta;
|
||||||
|
info = 43;
|
||||||
|
} // if
|
||||||
|
break;
|
||||||
|
|
||||||
|
case EstimatedParameter::measureErr_Corr:
|
||||||
|
k1 = estiParDesc.estParams[i].ID1;
|
||||||
|
k2 = estiParDesc.estParams[i].ID2;
|
||||||
|
// H(k1,k2) = xparam1(i)*sqrt(H(k1,k1)*H(k2,k2));
|
||||||
|
// H(k2,k1) = H(k1,k2);
|
||||||
|
H(k1, k2) = estParams(i)*sqrt(H(k1, k1)*H(k2, k2));
|
||||||
|
H(k2, k1) = H(k1, k2);
|
||||||
|
|
||||||
|
//[CholH,testH] = chol(H);
|
||||||
|
test = lapack::choleskyDecomp(H, "L");
|
||||||
|
if (test > 0)
|
||||||
|
{
|
||||||
|
mexPrintf("Caugth 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.
|
||||||
|
// We have to compute the eigenvalues of this matrix in order to build the penalty.
|
||||||
|
//a = diag(eig(H));
|
||||||
|
double delta = 0;
|
||||||
|
eigH.calculate(H); // get eigenvalues
|
||||||
|
//k = find(a < 0);
|
||||||
|
if (eigH.hasConverged())
|
||||||
|
{
|
||||||
|
const Vector &evH = eigH.getD();
|
||||||
|
for (i = 0; i < evH.getSize(); ++i)
|
||||||
|
if (evH(i) < 0)
|
||||||
|
delta -= evH(i);
|
||||||
|
}
|
||||||
|
logLikelihood = penalty+delta;
|
||||||
|
info = 44;
|
||||||
|
} // end if
|
||||||
|
break;
|
||||||
|
|
||||||
|
//if estim_params_.np > 0 // i.e. num of deep parameters >0
|
||||||
|
case EstimatedParameter::deepPar:
|
||||||
|
k = estiParDesc.estParams[i].ID1;
|
||||||
|
deepParams(k) = estParams(i);
|
||||||
|
break;
|
||||||
|
default:
|
||||||
|
logLikelihood = penalty;
|
||||||
|
info = 1;
|
||||||
|
} // end switch
|
||||||
|
} // end found
|
||||||
|
} //end for
|
||||||
|
};
|
||||||
|
|
||||||
|
|
||||||
};
|
};
|
||||||
|
|
||||||
|
|
|
@ -40,12 +40,6 @@ LogPosteriorDensity::LogPosteriorDensity(const std::string &modName, EstimatedPa
|
||||||
|
|
||||||
}
|
}
|
||||||
|
|
||||||
double
|
|
||||||
LogPosteriorDensity::compute(Matrix &steadyState, const Vector &estParams, Vector &deepParams, const MatrixConstView &data, Matrix &Q, Matrix &H, size_t presampleStart, int &info)
|
|
||||||
{
|
|
||||||
return -logLikelihoodMain.compute(steadyState, estParams, deepParams, data, Q, H, presampleStart, info)
|
|
||||||
-logPriorDensity.compute(estParams);
|
|
||||||
}
|
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* vector of log likelihoods for each Kalman step
|
* vector of log likelihoods for each Kalman step
|
||||||
|
|
|
@ -48,8 +48,15 @@ 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, int &info_arg);
|
double riccati_tol_arg, double lyapunov_tol_arg, int &info_arg);
|
||||||
double compute(Matrix &steadyState, const Vector &estParams, Vector &deepParams,
|
|
||||||
const MatrixConstView &data, Matrix &Q, Matrix &H, size_t presampleStart, int &info);
|
template <class VEC1, class VEC2>
|
||||||
|
double
|
||||||
|
compute(VEC1 &steadyState, VEC2 &estParams, VectorView &deepParams, const MatrixConstView &data, MatrixView &Q, Matrix &H, size_t presampleStart, int &info)
|
||||||
|
{
|
||||||
|
return -logLikelihoodMain.compute(steadyState, estParams, deepParams, data, Q, H, presampleStart, info)
|
||||||
|
-logPriorDensity.compute(estParams);
|
||||||
|
}
|
||||||
|
|
||||||
Vector&getLikVector();
|
Vector&getLikVector();
|
||||||
|
|
||||||
};
|
};
|
||||||
|
|
|
@ -33,20 +33,6 @@ LogPriorDensity::LogPriorDensity(EstimatedParametersDescription &estParsDesc_arg
|
||||||
{
|
{
|
||||||
};
|
};
|
||||||
|
|
||||||
double
|
|
||||||
LogPriorDensity::compute(const Vector &ep)
|
|
||||||
{
|
|
||||||
assert(estParsDesc.estParams.size() == ep.getSize());
|
|
||||||
double logPriorDensity = 0;
|
|
||||||
for (size_t i = 0; i < ep.getSize(); ++i)
|
|
||||||
{
|
|
||||||
logPriorDensity += log(((*(estParsDesc.estParams[i]).prior)).pdf(ep(i)));
|
|
||||||
if (std::isinf(fabs(logPriorDensity)))
|
|
||||||
return logPriorDensity;
|
|
||||||
}
|
|
||||||
return logPriorDensity;
|
|
||||||
};
|
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Return random number for prior fromits distribution
|
* Return random number for prior fromits distribution
|
||||||
*/
|
*/
|
||||||
|
|
|
@ -37,7 +37,20 @@ public:
|
||||||
LogPriorDensity(EstimatedParametersDescription &estParsDesc);
|
LogPriorDensity(EstimatedParametersDescription &estParsDesc);
|
||||||
virtual ~LogPriorDensity();
|
virtual ~LogPriorDensity();
|
||||||
|
|
||||||
double compute(const Vector &estParams);
|
template<class VEC>
|
||||||
|
double compute(VEC &ep)
|
||||||
|
{
|
||||||
|
assert(estParsDesc.estParams.size() == ep.getSize());
|
||||||
|
double logPriorDensity = 0;
|
||||||
|
for (size_t i = 0; i < ep.getSize(); ++i)
|
||||||
|
{
|
||||||
|
logPriorDensity += log(((*(estParsDesc.estParams[i]).prior)).pdf(ep(i)));
|
||||||
|
if (std::isinf(fabs(logPriorDensity)))
|
||||||
|
return logPriorDensity;
|
||||||
|
}
|
||||||
|
return logPriorDensity;
|
||||||
|
};
|
||||||
|
|
||||||
void computeNewParams(Vector &newParams);
|
void computeNewParams(Vector &newParams);
|
||||||
|
|
||||||
private:
|
private:
|
||||||
|
|
|
@ -51,41 +51,4 @@ ModelSolution::ModelSolution(const std::string &dynamicDllFile, size_t n_endo_a
|
||||||
back_inserter(zeta_back_mixed));
|
back_inserter(zeta_back_mixed));
|
||||||
}
|
}
|
||||||
|
|
||||||
void
|
|
||||||
ModelSolution::compute(VectorView &steadyState, const Vector &deepParams, Matrix &ghx, Matrix &ghu) throw (DecisionRules::BlanchardKahnException, GeneralizedSchurDecomposition::GSDException)
|
|
||||||
{
|
|
||||||
// compute Steady State
|
|
||||||
ComputeSteadyState(steadyState, deepParams);
|
|
||||||
|
|
||||||
// then get jacobian and
|
|
||||||
|
|
||||||
ComputeModelSolution(steadyState, deepParams, ghx, ghu);
|
|
||||||
|
|
||||||
}
|
|
||||||
|
|
||||||
void
|
|
||||||
ModelSolution::ComputeModelSolution(VectorView &steadyState, const Vector &deepParams, Matrix &ghx, Matrix &ghu) throw (DecisionRules::BlanchardKahnException, GeneralizedSchurDecomposition::GSDException)
|
|
||||||
{
|
|
||||||
// set extended Steady State
|
|
||||||
|
|
||||||
for (size_t i = 0; i < zeta_back_mixed.size(); i++)
|
|
||||||
llXsteadyState(i) = steadyState(zeta_back_mixed[i]);
|
|
||||||
|
|
||||||
for (size_t i = 0; i < n_endo; i++)
|
|
||||||
llXsteadyState(zeta_back_mixed.size() + i) = steadyState(i);
|
|
||||||
|
|
||||||
for (size_t i = 0; i < zeta_fwrd_mixed.size(); i++)
|
|
||||||
llXsteadyState(zeta_back_mixed.size() + n_endo + i) = steadyState(zeta_fwrd_mixed[i]);
|
|
||||||
|
|
||||||
//get jacobian
|
|
||||||
dynamicDLLp.eval(llXsteadyState, Mx, deepParams, steadyState, residual, &jacobian, NULL, NULL);
|
|
||||||
|
|
||||||
//compute rules
|
|
||||||
decisionRules.compute(jacobian, ghx, ghu);
|
|
||||||
}
|
|
||||||
void
|
|
||||||
ModelSolution::ComputeSteadyState(VectorView &steadyState, const Vector &deepParams)
|
|
||||||
{
|
|
||||||
// does nothig for time being.
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
|
@ -42,7 +42,17 @@ public:
|
||||||
const std::vector<size_t> &zeta_back_arg, const std::vector<size_t> &zeta_mixed_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> &zeta_static_arg, double qz_criterium);
|
||||||
virtual ~ModelSolution() {};
|
virtual ~ModelSolution() {};
|
||||||
void compute(VectorView &steadyState, const Vector &deepParams, Matrix &ghx, Matrix &ghu) throw (DecisionRules::BlanchardKahnException, GeneralizedSchurDecomposition::GSDException);
|
template <class VEC>
|
||||||
|
void compute(VEC &steadyState, const VectorView &deepParams, Matrix &ghx, Matrix &ghu) throw (DecisionRules::BlanchardKahnException, GeneralizedSchurDecomposition::GSDException)
|
||||||
|
{
|
||||||
|
// compute Steady State
|
||||||
|
ComputeSteadyState(steadyState, deepParams);
|
||||||
|
|
||||||
|
// then get jacobian and
|
||||||
|
|
||||||
|
ComputeModelSolution(steadyState, deepParams, ghx, ghu);
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
private:
|
private:
|
||||||
const size_t n_endo;
|
const size_t n_endo;
|
||||||
|
@ -56,8 +66,34 @@ private:
|
||||||
DynamicModelDLL dynamicDLLp;
|
DynamicModelDLL dynamicDLLp;
|
||||||
Vector llXsteadyState;
|
Vector llXsteadyState;
|
||||||
//Matrix jacobian;
|
//Matrix jacobian;
|
||||||
void ComputeModelSolution(VectorView &steadyState, const Vector &deepParams, Matrix &ghx, Matrix &ghu) throw (DecisionRules::BlanchardKahnException, GeneralizedSchurDecomposition::GSDException);
|
template <class VEC>
|
||||||
void ComputeSteadyState(VectorView &steadyState, const Vector &deepParams);
|
void ComputeModelSolution(VEC &steadyState, const VectorView &deepParams,
|
||||||
|
Matrix &ghx, Matrix &ghu)
|
||||||
|
throw (DecisionRules::BlanchardKahnException, GeneralizedSchurDecomposition::GSDException)
|
||||||
|
{
|
||||||
|
// set extended Steady State
|
||||||
|
|
||||||
|
for (size_t i = 0; i < zeta_back_mixed.size(); i++)
|
||||||
|
llXsteadyState(i) = steadyState(zeta_back_mixed[i]);
|
||||||
|
|
||||||
|
for (size_t i = 0; i < n_endo; i++)
|
||||||
|
llXsteadyState(zeta_back_mixed.size() + i) = steadyState(i);
|
||||||
|
|
||||||
|
for (size_t i = 0; i < zeta_fwrd_mixed.size(); i++)
|
||||||
|
llXsteadyState(zeta_back_mixed.size() + n_endo + i) = steadyState(zeta_fwrd_mixed[i]);
|
||||||
|
|
||||||
|
//get jacobian
|
||||||
|
dynamicDLLp.eval(llXsteadyState, Mx, deepParams, steadyState, residual, &jacobian, NULL, NULL);
|
||||||
|
|
||||||
|
//compute rules
|
||||||
|
decisionRules.compute(jacobian, ghx, ghu);
|
||||||
|
}
|
||||||
|
template <class VEC>
|
||||||
|
void ComputeSteadyState(VEC &steadyState, const VectorView &deepParams)
|
||||||
|
{
|
||||||
|
// does nothig for time being.
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
};
|
};
|
||||||
|
|
||||||
|
|
|
@ -1,92 +0,0 @@
|
||||||
/*
|
|
||||||
* Copyright (C) 2010-2011 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/>.
|
|
||||||
*/
|
|
||||||
|
|
||||||
#include "RandomWalkMetropolisHastings.hh"
|
|
||||||
|
|
||||||
#include <iostream>
|
|
||||||
#include <fstream>
|
|
||||||
|
|
||||||
double
|
|
||||||
RandomWalkMetropolisHastings::compute(VectorView &mhLogPostDens, MatrixView &mhParams, Matrix &steadyState,
|
|
||||||
Vector &estParams, Vector &deepParams, const MatrixConstView &data, Matrix &Q, Matrix &H,
|
|
||||||
const size_t presampleStart, int &info, const size_t startDraw, size_t nMHruns,
|
|
||||||
LogPosteriorDensity &lpd, Proposal &pDD, EstimatedParametersDescription &epd)
|
|
||||||
{
|
|
||||||
//streambuf *likbuf, *drawbuf *backup;
|
|
||||||
std::ofstream urandfilestr, drawfilestr;
|
|
||||||
urandfilestr.open("urand.csv");
|
|
||||||
drawfilestr.open("paramdraws.csv");
|
|
||||||
|
|
||||||
bool overbound;
|
|
||||||
double newLogpost, logpost, urand;
|
|
||||||
size_t count, accepted = 0;
|
|
||||||
parDraw = estParams;
|
|
||||||
|
|
||||||
logpost = -lpd.compute(steadyState, estParams, deepParams, data, Q, H, presampleStart, info);
|
|
||||||
|
|
||||||
for (size_t run = startDraw - 1; run < nMHruns; ++run)
|
|
||||||
{
|
|
||||||
overbound = false;
|
|
||||||
pDD.draw(parDraw, newParDraw);
|
|
||||||
for (count = 0; count < parDraw.getSize(); ++count)
|
|
||||||
{
|
|
||||||
overbound = (newParDraw(count) < epd.estParams[count].lower_bound || newParDraw(count) > epd.estParams[count].upper_bound);
|
|
||||||
if (overbound)
|
|
||||||
{
|
|
||||||
newLogpost = -INFINITY;
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
if (!overbound)
|
|
||||||
{
|
|
||||||
try
|
|
||||||
{
|
|
||||||
newLogpost = -lpd.compute(steadyState, newParDraw, deepParams, data, Q, H, presampleStart, info);
|
|
||||||
}
|
|
||||||
catch (const std::exception &e)
|
|
||||||
{
|
|
||||||
throw; // for now handle the system and other errors higher-up
|
|
||||||
}
|
|
||||||
catch (...)
|
|
||||||
{
|
|
||||||
newLogpost = -INFINITY;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
urand = pDD.selectionTestDraw();
|
|
||||||
if ((newLogpost > -INFINITY) && log(urand) < newLogpost-logpost)
|
|
||||||
{
|
|
||||||
parDraw = newParDraw;
|
|
||||||
logpost = newLogpost;
|
|
||||||
accepted++;
|
|
||||||
}
|
|
||||||
mat::get_row(mhParams, run) = parDraw;
|
|
||||||
mhLogPostDens(run) = logpost;
|
|
||||||
|
|
||||||
urandfilestr << urand << "\n"; //","
|
|
||||||
for (size_t c = 0; c < newParDraw.getSize()-1; ++c)
|
|
||||||
drawfilestr << newParDraw(c) << ",";
|
|
||||||
drawfilestr << newParDraw(newParDraw.getSize()-1) << "\n";
|
|
||||||
}
|
|
||||||
|
|
||||||
urandfilestr.close();
|
|
||||||
drawfilestr.close();
|
|
||||||
|
|
||||||
return (double) accepted/(nMHruns-startDraw+1);
|
|
||||||
}
|
|
||||||
|
|
|
@ -26,6 +26,8 @@
|
||||||
#if !defined(A6BBC5E0_598E_4863_B7FF_E87320056B80__INCLUDED_)
|
#if !defined(A6BBC5E0_598E_4863_B7FF_E87320056B80__INCLUDED_)
|
||||||
#define A6BBC5E0_598E_4863_B7FF_E87320056B80__INCLUDED_
|
#define A6BBC5E0_598E_4863_B7FF_E87320056B80__INCLUDED_
|
||||||
|
|
||||||
|
#include <iostream>
|
||||||
|
#include <fstream>
|
||||||
#include "LogPosteriorDensity.hh"
|
#include "LogPosteriorDensity.hh"
|
||||||
#include "Proposal.hh"
|
#include "Proposal.hh"
|
||||||
|
|
||||||
|
@ -41,11 +43,75 @@ public:
|
||||||
{
|
{
|
||||||
};
|
};
|
||||||
virtual ~RandomWalkMetropolisHastings() {};
|
virtual ~RandomWalkMetropolisHastings() {};
|
||||||
virtual double compute(VectorView &mhLogPostDens, MatrixView &mhParams, Matrix &steadyState,
|
|
||||||
Vector &estParams, Vector &deepParams, const MatrixConstView &data, Matrix &Q, Matrix &H,
|
template<class VEC1>
|
||||||
const size_t presampleStart, int &info, const size_t startDraw, size_t nMHruns,
|
double compute(VectorView &mhLogPostDens, MatrixView &mhParams, VEC1 &steadyState,
|
||||||
LogPosteriorDensity &logPosteriorDensity, Proposal &proposalDrawDistribution,
|
Vector &estParams, VectorView &deepParams, const MatrixConstView &data, MatrixView &Q, Matrix &H,
|
||||||
EstimatedParametersDescription &epd);
|
const size_t presampleStart, int &info, const size_t startDraw, size_t nMHruns,
|
||||||
|
LogPosteriorDensity &lpd, Proposal &pDD, EstimatedParametersDescription &epd)
|
||||||
|
{
|
||||||
|
//streambuf *likbuf, *drawbuf *backup;
|
||||||
|
std::ofstream urandfilestr, drawfilestr;
|
||||||
|
urandfilestr.open("urand.csv");
|
||||||
|
drawfilestr.open("paramdraws.csv");
|
||||||
|
|
||||||
|
bool overbound;
|
||||||
|
double newLogpost, logpost, urand;
|
||||||
|
size_t count, accepted = 0;
|
||||||
|
parDraw = estParams;
|
||||||
|
|
||||||
|
logpost = -lpd.compute(steadyState, estParams, deepParams, data, Q, H, presampleStart, info);
|
||||||
|
|
||||||
|
for (size_t run = startDraw - 1; run < nMHruns; ++run)
|
||||||
|
{
|
||||||
|
overbound = false;
|
||||||
|
pDD.draw(parDraw, newParDraw);
|
||||||
|
for (count = 0; count < parDraw.getSize(); ++count)
|
||||||
|
{
|
||||||
|
overbound = (newParDraw(count) < epd.estParams[count].lower_bound || newParDraw(count) > epd.estParams[count].upper_bound);
|
||||||
|
if (overbound)
|
||||||
|
{
|
||||||
|
newLogpost = -INFINITY;
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
if (!overbound)
|
||||||
|
{
|
||||||
|
try
|
||||||
|
{
|
||||||
|
newLogpost = -lpd.compute(steadyState, newParDraw, deepParams, data, Q, H, presampleStart, info);
|
||||||
|
}
|
||||||
|
catch (const std::exception &e)
|
||||||
|
{
|
||||||
|
throw; // for now handle the system and other errors higher-up
|
||||||
|
}
|
||||||
|
catch (...)
|
||||||
|
{
|
||||||
|
newLogpost = -INFINITY;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
urand = pDD.selectionTestDraw();
|
||||||
|
if ((newLogpost > -INFINITY) && log(urand) < newLogpost-logpost)
|
||||||
|
{
|
||||||
|
parDraw = newParDraw;
|
||||||
|
logpost = newLogpost;
|
||||||
|
accepted++;
|
||||||
|
}
|
||||||
|
mat::get_row(mhParams, run) = parDraw;
|
||||||
|
mhLogPostDens(run) = logpost;
|
||||||
|
|
||||||
|
urandfilestr << urand << "\n"; //","
|
||||||
|
for (size_t c = 0; c < newParDraw.getSize()-1; ++c)
|
||||||
|
drawfilestr << newParDraw(c) << ",";
|
||||||
|
drawfilestr << newParDraw(newParDraw.getSize()-1) << "\n";
|
||||||
|
}
|
||||||
|
|
||||||
|
urandfilestr.close();
|
||||||
|
drawfilestr.close();
|
||||||
|
|
||||||
|
return (double) accepted/(nMHruns-startDraw+1);
|
||||||
|
};
|
||||||
|
|
||||||
};
|
};
|
||||||
|
|
||||||
#endif // !defined(A6BBC5E0_598E_4863_B7FF_E87320056B80__INCLUDED_)
|
#endif // !defined(A6BBC5E0_598E_4863_B7FF_E87320056B80__INCLUDED_)
|
||||||
|
|
|
@ -114,8 +114,8 @@ fillEstParamsInfo(const mxArray *estim_params_info, EstimatedParameter::pType ty
|
||||||
|
|
||||||
int
|
int
|
||||||
sampleMHMC(LogPosteriorDensity &lpd, RandomWalkMetropolisHastings &rwmh,
|
sampleMHMC(LogPosteriorDensity &lpd, RandomWalkMetropolisHastings &rwmh,
|
||||||
Matrix &steadyState, Vector &estParams, Vector &deepParams, const MatrixConstView &data,
|
VectorView &steadyState, VectorConstView &estParams, VectorView &deepParams, const MatrixConstView &data,
|
||||||
Matrix &Q, Matrix &H, size_t presampleStart, int &info, const VectorConstView &nruns,
|
MatrixView &Q, Matrix &H, size_t presampleStart, int &info, 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)
|
||||||
{
|
{
|
||||||
|
@ -539,26 +539,26 @@ sampleMHMC(LogPosteriorDensity &lpd, RandomWalkMetropolisHastings &rwmh,
|
||||||
mhFName = ssFName.str();
|
mhFName = ssFName.str();
|
||||||
fidlog = fopen(mhFName.c_str(), "a");
|
fidlog = fopen(mhFName.c_str(), "a");
|
||||||
fprintf(fidlog, "\n");
|
fprintf(fidlog, "\n");
|
||||||
fprintf(fidlog, "%% Mh%dBlck%d ( %s %s )\n", (int) NewFileVw(b-1), b, __DATE__, __TIME__);
|
fprintf(fidlog, "%% Mh%dBlck%lu ( %s %s )\n", (int) NewFileVw(b-1), b, __DATE__, __TIME__);
|
||||||
fprintf(fidlog, " \n");
|
fprintf(fidlog, " \n");
|
||||||
fprintf(fidlog, " Number of simulations.: %d \n", currInitSizeArray); // (length(logpo2)) ');
|
fprintf(fidlog, " Number of simulations.: %lu \n", currInitSizeArray); // (length(logpo2)) ');
|
||||||
fprintf(fidlog, " Acceptation rate......: %f \n", jsux);
|
fprintf(fidlog, " Acceptation rate......: %f \n", jsux);
|
||||||
fprintf(fidlog, " Posterior mean........:\n");
|
fprintf(fidlog, " Posterior mean........:\n");
|
||||||
for (size_t i = 0; i < npar; ++i)
|
for (size_t i = 0; i < npar; ++i)
|
||||||
{
|
{
|
||||||
VectorView mhpdColVw = mat::get_col(mhParamDraws, i);
|
VectorView mhpdColVw = mat::get_col(mhParamDraws, i);
|
||||||
fprintf(fidlog, " params: %d : %f \n", i+1, vec::meanSumMinMax(dsum, dmin, dmax, mhpdColVw));
|
fprintf(fidlog, " params: %lu : %f \n", i+1, vec::meanSumMinMax(dsum, dmin, dmax, mhpdColVw));
|
||||||
MinMax(i, iMin) = dmin;
|
MinMax(i, iMin) = dmin;
|
||||||
MinMax(i, iMax) = dmax;
|
MinMax(i, iMax) = dmax;
|
||||||
} // end
|
} // end
|
||||||
fprintf(fidlog, " log2po: %f \n", vec::meanSumMinMax(dsum, dmin, dmax, mhLogPostDens));
|
fprintf(fidlog, " log2po: %f \n", vec::meanSumMinMax(dsum, dmin, dmax, mhLogPostDens));
|
||||||
fprintf(fidlog, " Minimum value.........:\n");;
|
fprintf(fidlog, " Minimum value.........:\n");;
|
||||||
for (size_t i = 0; i < npar; ++i)
|
for (size_t i = 0; i < npar; ++i)
|
||||||
fprintf(fidlog, " params: %d : %f \n", i+1, MinMax(i, iMin));
|
fprintf(fidlog, " params: %lu : %f \n", i+1, MinMax(i, iMin));
|
||||||
fprintf(fidlog, " log2po: %f \n", dmin);
|
fprintf(fidlog, " log2po: %f \n", dmin);
|
||||||
fprintf(fidlog, " Maximum value.........:\n");
|
fprintf(fidlog, " Maximum value.........:\n");
|
||||||
for (size_t i = 0; i < npar; ++i)
|
for (size_t i = 0; i < npar; ++i)
|
||||||
fprintf(fidlog, " params: %d : %f \n", i+1, MinMax(i, iMax));
|
fprintf(fidlog, " params: %lu : %f \n", i+1, MinMax(i, iMax));
|
||||||
fprintf(fidlog, " log2po: %f \n", dmax);
|
fprintf(fidlog, " log2po: %f \n", dmax);
|
||||||
fprintf(fidlog, " \n");
|
fprintf(fidlog, " \n");
|
||||||
fclose(fidlog);
|
fclose(fidlog);
|
||||||
|
@ -624,12 +624,12 @@ sampleMHMC(LogPosteriorDensity &lpd, RandomWalkMetropolisHastings &rwmh,
|
||||||
}
|
}
|
||||||
|
|
||||||
int
|
int
|
||||||
logMCMCposterior(const VectorConstView &estParams, const MatrixConstView &data,
|
logMCMCposterior(VectorConstView &estParams, const MatrixConstView &data,
|
||||||
const size_t fblock, const size_t nBlocks, const VectorConstView &nMHruns, const MatrixConstView &D)
|
const size_t fblock, const size_t nBlocks, const VectorConstView &nMHruns, const MatrixConstView &D,
|
||||||
|
VectorView &steadyState, VectorView &deepParams, MatrixView &Q, Matrix &H)
|
||||||
{
|
{
|
||||||
// Retrieve pointers to global variables
|
// Retrieve pointers to global variables
|
||||||
const mxArray *M_ = mexGetVariablePtr("global", "M_");
|
const mxArray *M_ = mexGetVariablePtr("global", "M_");
|
||||||
const mxArray *oo_ = mexGetVariablePtr("global", "oo_");
|
|
||||||
const mxArray *options_ = mexGetVariablePtr("global", "options_");
|
const mxArray *options_ = mexGetVariablePtr("global", "options_");
|
||||||
const mxArray *estim_params_ = mexGetVariablePtr("global", "estim_params_");
|
const mxArray *estim_params_ = mexGetVariablePtr("global", "estim_params_");
|
||||||
|
|
||||||
|
@ -643,7 +643,7 @@ logMCMCposterior(const VectorConstView &estParams, const MatrixConstView &data,
|
||||||
|
|
||||||
size_t n_endo = (size_t) *mxGetPr(mxGetField(M_, 0, "endo_nbr"));
|
size_t n_endo = (size_t) *mxGetPr(mxGetField(M_, 0, "endo_nbr"));
|
||||||
size_t n_exo = (size_t) *mxGetPr(mxGetField(M_, 0, "exo_nbr"));
|
size_t n_exo = (size_t) *mxGetPr(mxGetField(M_, 0, "exo_nbr"));
|
||||||
size_t n_param = (size_t) *mxGetPr(mxGetField(M_, 0, "param_nbr"));
|
|
||||||
size_t n_estParams = estParams.getSize();
|
size_t n_estParams = estParams.getSize();
|
||||||
|
|
||||||
std::vector<size_t> zeta_fwrd, zeta_back, zeta_mixed, zeta_static;
|
std::vector<size_t> zeta_fwrd, zeta_back, zeta_mixed, zeta_static;
|
||||||
|
@ -702,26 +702,8 @@ logMCMCposterior(const VectorConstView &estParams, const MatrixConstView &data,
|
||||||
LogPosteriorDensity lpd(dynamicDllFile, epd, n_endo, n_exo, zeta_fwrd, zeta_back, zeta_mixed, zeta_static,
|
LogPosteriorDensity lpd(dynamicDllFile, epd, n_endo, n_exo, zeta_fwrd, zeta_back, zeta_mixed, zeta_static,
|
||||||
qz_criterium, varobs, riccati_tol, lyapunov_tol, info);
|
qz_criterium, varobs, riccati_tol, lyapunov_tol, info);
|
||||||
|
|
||||||
// Construct arguments of compute() method
|
|
||||||
Matrix steadyState(n_endo, 1);
|
|
||||||
mat::get_col(steadyState, 0) = VectorConstView(mxGetPr(mxGetField(oo_, 0, "steady_state")), n_endo, 1);
|
|
||||||
|
|
||||||
Vector estParams2(n_estParams);
|
|
||||||
estParams2 = estParams;
|
|
||||||
Vector deepParams(n_param);
|
|
||||||
deepParams = VectorConstView(mxGetPr(mxGetField(M_, 0, "params")), n_param, 1);
|
|
||||||
Matrix Q(n_exo);
|
|
||||||
Q = MatrixConstView(mxGetPr(mxGetField(M_, 0, "Sigma_e")), n_exo, n_exo, n_exo);
|
|
||||||
|
|
||||||
Matrix H(n_varobs);
|
|
||||||
const mxArray *H_mx = mxGetField(M_, 0, "H");
|
|
||||||
if (mxGetM(H_mx) == 1 && mxGetN(H_mx) == 1 && *mxGetPr(H_mx) == 0)
|
|
||||||
H.setAll(0.0);
|
|
||||||
else
|
|
||||||
H = MatrixConstView(mxGetPr(mxGetField(M_, 0, "H")), n_varobs, n_varobs, n_varobs);
|
|
||||||
|
|
||||||
// Construct MHMCMC Sampler
|
// Construct MHMCMC Sampler
|
||||||
RandomWalkMetropolisHastings rwmh(estParams2.getSize());
|
RandomWalkMetropolisHastings rwmh(estParams.getSize());
|
||||||
// Construct GaussianPrior drawDistribution m=0, sd=1
|
// Construct GaussianPrior drawDistribution m=0, sd=1
|
||||||
GaussianPrior drawGaussDist01(0.0, 1.0, -INFINITY, INFINITY, 0.0, 1.0);
|
GaussianPrior drawGaussDist01(0.0, 1.0, -INFINITY, INFINITY, 0.0, 1.0);
|
||||||
// get Jscale = diag(bayestopt_.jscale);
|
// get Jscale = diag(bayestopt_.jscale);
|
||||||
|
@ -731,7 +713,7 @@ logMCMCposterior(const 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, estParams2, deepParams, data, Q, H, presample, info,
|
int lastMHblockArrayLine = sampleMHMC(lpd, rwmh, steadyState, estParams, deepParams, data, Q, H, presample, info,
|
||||||
nMHruns, fblock, nBlocks, pdd, epd, resultsFileStem, console_mode, load_mh_file);
|
nMHruns, fblock, nBlocks, pdd, epd, resultsFileStem, console_mode, load_mh_file);
|
||||||
|
|
||||||
// Cleanups
|
// Cleanups
|
||||||
|
@ -746,8 +728,8 @@ void
|
||||||
mexFunction(int nlhs, mxArray *plhs[],
|
mexFunction(int nlhs, mxArray *plhs[],
|
||||||
int nrhs, const mxArray *prhs[])
|
int nrhs, const mxArray *prhs[])
|
||||||
{
|
{
|
||||||
if (nrhs != 6)
|
if (nrhs != 10)
|
||||||
DYN_MEX_FUNC_ERR_MSG_TXT("logposterior: exactly six arguments are required.");
|
DYN_MEX_FUNC_ERR_MSG_TXT("logposterior: exactly 11 arguments are required.");
|
||||||
if (nlhs != 2)
|
if (nlhs != 2)
|
||||||
DYN_MEX_FUNC_ERR_MSG_TXT("logposterior: exactly two return arguments are required.");
|
DYN_MEX_FUNC_ERR_MSG_TXT("logposterior: exactly two return arguments are required.");
|
||||||
|
|
||||||
|
@ -761,18 +743,42 @@ mexFunction(int nlhs, mxArray *plhs[],
|
||||||
if (!mxIsDouble(prhs[1]))
|
if (!mxIsDouble(prhs[1]))
|
||||||
DYN_MEX_FUNC_ERR_MSG_TXT("logposterior: Second argument must be a matrix of double-precision numbers");
|
DYN_MEX_FUNC_ERR_MSG_TXT("logposterior: Second argument must be a matrix of double-precision numbers");
|
||||||
|
|
||||||
MatrixConstView data(mxGetPr(prhs[1]), mxGetM(prhs[1]), mxGetN(prhs[1]), mxGetM(prhs[1]));
|
|
||||||
|
|
||||||
size_t fblock = (size_t) mxGetScalar(prhs[2]);
|
size_t fblock = (size_t) mxGetScalar(prhs[2]);
|
||||||
size_t nBlocks = (size_t) mxGetScalar(prhs[3]);
|
size_t nBlocks = (size_t) mxGetScalar(prhs[3]);
|
||||||
VectorConstView nMHruns(mxGetPr(prhs[4]), mxGetM(prhs[4]), 1);
|
VectorConstView nMHruns(mxGetPr(prhs[4]), mxGetM(prhs[4]), 1);
|
||||||
|
MatrixConstView D(mxGetPr(prhs[5]), mxGetM(prhs[5]), mxGetN(prhs[5]), mxGetM(prhs[5]));
|
||||||
|
const mxArray *dataset = prhs[6];
|
||||||
|
const mxArray *options_ = prhs[7];
|
||||||
|
const mxArray *M_ = prhs[8];
|
||||||
|
const mxArray *bayestopt_ = prhs[9];
|
||||||
|
const mxArray *oo_ = prhs[10];
|
||||||
|
|
||||||
assert(nMHruns.getSize() == nBlocks);
|
assert(nMHruns.getSize() == nBlocks);
|
||||||
|
|
||||||
MatrixConstView D(mxGetPr(prhs[5]), mxGetM(prhs[5]), mxGetN(prhs[5]), mxGetM(prhs[5]));
|
mxArray *dataset_data = mxGetField(dataset,0,"data");
|
||||||
|
MatrixConstView data(mxGetPr(dataset_data), mxGetM(dataset_data), mxGetN(dataset_data), mxGetM(dataset_data));
|
||||||
|
|
||||||
|
int endo_nbr = *(int*)mxGetPr(mxGetField(M_, 0, "endo_nbr"));
|
||||||
|
int exo_nbr = *(int*)mxGetPr(mxGetField(M_, 0, "exo_nbr"));
|
||||||
|
int param_nbr = *(int*)mxGetPr(mxGetField(M_, 0, "param_nbr"));
|
||||||
|
int varobs_nbr = mxGetN(mxGetField(options_, 0, "varobs"));
|
||||||
|
|
||||||
|
VectorView steadyState(mxGetPr(mxGetField(oo_,0,"steady_state")),endo_nbr, 1);
|
||||||
|
VectorView deepParams(mxGetPr(mxGetField(M_, 0, "params")),param_nbr,1);
|
||||||
|
|
||||||
|
MatrixView Q(mxGetPr(mxGetField(M_, 0, "Sigma_e")), exo_nbr, exo_nbr, exo_nbr);
|
||||||
|
|
||||||
|
Matrix H(varobs_nbr,varobs_nbr);
|
||||||
|
const mxArray *H_mx = mxGetField(M_, 0, "H");
|
||||||
|
if (mxGetM(H_mx) == 1 && mxGetN(H_mx) == 1 && *mxGetPr(H_mx) == 0)
|
||||||
|
H.setAll(0.0);
|
||||||
|
else
|
||||||
|
H = MatrixConstView(mxGetPr(H_mx), varobs_nbr, varobs_nbr, varobs_nbr);
|
||||||
|
|
||||||
//calculate MHMCMC draws and get get last line run in the last MH block sub-array
|
//calculate MHMCMC draws and get get last line run in the last MH block sub-array
|
||||||
try
|
try
|
||||||
{
|
{
|
||||||
int lastMHblockArrayLine = logMCMCposterior(estParams, data, fblock, nBlocks, nMHruns, D);
|
int lastMHblockArrayLine = logMCMCposterior(estParams, data, fblock, nBlocks, nMHruns, D, steadyState, deepParams, Q, H);
|
||||||
plhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL);
|
plhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL);
|
||||||
*mxGetPr(plhs[1]) = (double) lastMHblockArrayLine;
|
*mxGetPr(plhs[1]) = (double) lastMHblockArrayLine;
|
||||||
}
|
}
|
||||||
|
|
|
@ -21,6 +21,7 @@
|
||||||
#include <vector>
|
#include <vector>
|
||||||
#include <algorithm>
|
#include <algorithm>
|
||||||
#include <functional>
|
#include <functional>
|
||||||
|
#include <sstream>
|
||||||
|
|
||||||
#include "Vector.hh"
|
#include "Vector.hh"
|
||||||
#include "Matrix.hh"
|
#include "Matrix.hh"
|
||||||
|
@ -39,11 +40,10 @@ public:
|
||||||
};
|
};
|
||||||
|
|
||||||
void
|
void
|
||||||
fillEstParamsInfo(const mxArray *estim_params_info, EstimatedParameter::pType type,
|
fillEstParamsInfo(const mxArray *bayestopt_, const mxArray *estim_params_info, EstimatedParameter::pType type,
|
||||||
std::vector<EstimatedParameter> &estParamsInfo)
|
std::vector<EstimatedParameter> &estParamsInfo)
|
||||||
{
|
{
|
||||||
// execute once only
|
// execute once only
|
||||||
static const mxArray *bayestopt_ = mexGetVariablePtr("global", "bayestopt_");
|
|
||||||
static const mxArray *bayestopt_ubp = mxGetField(bayestopt_, 0, "ub"); // upper bound
|
static const mxArray *bayestopt_ubp = mxGetField(bayestopt_, 0, "ub"); // upper bound
|
||||||
static const mxArray *bayestopt_lbp = mxGetField(bayestopt_, 0, "lb"); // lower bound
|
static const mxArray *bayestopt_lbp = mxGetField(bayestopt_, 0, "lb"); // lower bound
|
||||||
static const mxArray *bayestopt_p1p = mxGetField(bayestopt_, 0, "p1"); // prior mean
|
static const mxArray *bayestopt_p1p = mxGetField(bayestopt_, 0, "p1"); // prior mean
|
||||||
|
@ -100,15 +100,13 @@ fillEstParamsInfo(const mxArray *estim_params_info, EstimatedParameter::pType ty
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
template <class VEC1, class VEC2>
|
||||||
double
|
double
|
||||||
logposterior(const VectorConstView &estParams, const MatrixConstView &data)
|
logposterior(VEC1 &estParams, const MatrixConstView &data,
|
||||||
|
const mxArray *options_, const mxArray *M_, const mxArray *estim_params_,
|
||||||
|
const mxArray *bayestopt_, const mxArray *oo_, VEC2 &steadyState, double *trend_coeff,
|
||||||
|
int &info, VectorView &deepParams, Matrix &H, MatrixView &Q)
|
||||||
{
|
{
|
||||||
// Retrieve pointers to global variables
|
|
||||||
const mxArray *M_ = mexGetVariablePtr("global", "M_");
|
|
||||||
const mxArray *oo_ = mexGetVariablePtr("global", "oo_");
|
|
||||||
const mxArray *options_ = mexGetVariablePtr("global", "options_");
|
|
||||||
const mxArray *estim_params_ = mexGetVariablePtr("global", "estim_params_");
|
|
||||||
|
|
||||||
// Construct arguments of constructor of LogLikelihoodMain
|
// Construct arguments of constructor of LogLikelihoodMain
|
||||||
char *fName = mxArrayToString(mxGetField(M_, 0, "fname"));
|
char *fName = mxArrayToString(mxGetField(M_, 0, "fname"));
|
||||||
std::string dynamicDllFile(fName);
|
std::string dynamicDllFile(fName);
|
||||||
|
@ -150,6 +148,7 @@ logposterior(const VectorConstView &estParams, const MatrixConstView &data)
|
||||||
throw LogposteriorMexErrMsgTxtException("options_.varobs_id must be a row vector");
|
throw LogposteriorMexErrMsgTxtException("options_.varobs_id must be a row vector");
|
||||||
|
|
||||||
size_t n_varobs = mxGetN(varobs_mx);
|
size_t n_varobs = mxGetN(varobs_mx);
|
||||||
|
// substract 1.0 from obsverved variables index
|
||||||
std::transform(mxGetPr(varobs_mx), mxGetPr(varobs_mx) + n_varobs, back_inserter(varobs),
|
std::transform(mxGetPr(varobs_mx), mxGetPr(varobs_mx) + n_varobs, back_inserter(varobs),
|
||||||
std::bind2nd(std::minus<size_t>(), 1));
|
std::bind2nd(std::minus<size_t>(), 1));
|
||||||
|
|
||||||
|
@ -160,44 +159,27 @@ logposterior(const VectorConstView &estParams, const MatrixConstView &data)
|
||||||
estSubsamples.push_back(EstimationSubsample(0, data.getCols() - 1));
|
estSubsamples.push_back(EstimationSubsample(0, data.getCols() - 1));
|
||||||
|
|
||||||
std::vector<EstimatedParameter> estParamsInfo;
|
std::vector<EstimatedParameter> estParamsInfo;
|
||||||
fillEstParamsInfo(mxGetField(estim_params_, 0, "var_exo"), EstimatedParameter::shock_SD,
|
fillEstParamsInfo(bayestopt_, mxGetField(estim_params_, 0, "var_exo"), EstimatedParameter::shock_SD,
|
||||||
estParamsInfo);
|
estParamsInfo);
|
||||||
fillEstParamsInfo(mxGetField(estim_params_, 0, "var_endo"), EstimatedParameter::measureErr_SD,
|
fillEstParamsInfo(bayestopt_, mxGetField(estim_params_, 0, "var_endo"), EstimatedParameter::measureErr_SD,
|
||||||
estParamsInfo);
|
estParamsInfo);
|
||||||
fillEstParamsInfo(mxGetField(estim_params_, 0, "corrx"), EstimatedParameter::shock_Corr,
|
fillEstParamsInfo(bayestopt_, mxGetField(estim_params_, 0, "corrx"), EstimatedParameter::shock_Corr,
|
||||||
estParamsInfo);
|
estParamsInfo);
|
||||||
fillEstParamsInfo(mxGetField(estim_params_, 0, "corrn"), EstimatedParameter::measureErr_Corr,
|
fillEstParamsInfo(bayestopt_, mxGetField(estim_params_, 0, "corrn"), EstimatedParameter::measureErr_Corr,
|
||||||
estParamsInfo);
|
estParamsInfo);
|
||||||
fillEstParamsInfo(mxGetField(estim_params_, 0, "param_vals"), EstimatedParameter::deepPar,
|
fillEstParamsInfo(bayestopt_, mxGetField(estim_params_, 0, "param_vals"), EstimatedParameter::deepPar,
|
||||||
estParamsInfo);
|
estParamsInfo);
|
||||||
|
|
||||||
EstimatedParametersDescription epd(estSubsamples, estParamsInfo);
|
EstimatedParametersDescription epd(estSubsamples, estParamsInfo);
|
||||||
|
|
||||||
// Allocate LogPosteriorDensity object
|
// Allocate LogPosteriorDensity object
|
||||||
int info;
|
|
||||||
LogPosteriorDensity lpd(dynamicDllFile, epd, n_endo, n_exo, zeta_fwrd, zeta_back, zeta_mixed, zeta_static,
|
LogPosteriorDensity lpd(dynamicDllFile, epd, n_endo, n_exo, zeta_fwrd, zeta_back, zeta_mixed, zeta_static,
|
||||||
qz_criterium, varobs, riccati_tol, lyapunov_tol, info);
|
qz_criterium, varobs, riccati_tol, lyapunov_tol, info);
|
||||||
|
|
||||||
// Construct arguments of compute() method
|
// Construct arguments of compute() method
|
||||||
Matrix steadyState(n_endo, 1);
|
|
||||||
mat::get_col(steadyState, 0) = VectorConstView(mxGetPr(mxGetField(oo_, 0, "steady_state")), n_endo, 1);
|
|
||||||
|
|
||||||
Vector estParams2(n_estParams);
|
|
||||||
estParams2 = estParams;
|
|
||||||
Vector deepParams(n_param);
|
|
||||||
deepParams = VectorConstView(mxGetPr(mxGetField(M_, 0, "params")), n_param, 1);
|
|
||||||
Matrix Q(n_exo);
|
|
||||||
Q = MatrixConstView(mxGetPr(mxGetField(M_, 0, "Sigma_e")), n_exo, n_exo, n_exo);
|
|
||||||
|
|
||||||
Matrix H(n_varobs);
|
|
||||||
const mxArray *H_mx = mxGetField(M_, 0, "H");
|
|
||||||
if (mxGetM(H_mx) == 1 && mxGetN(H_mx) == 1 && *mxGetPr(H_mx) == 0)
|
|
||||||
H.setAll(0.0);
|
|
||||||
else
|
|
||||||
H = MatrixConstView(mxGetPr(mxGetField(M_, 0, "H")), n_varobs, n_varobs, n_varobs);
|
|
||||||
|
|
||||||
// Compute the posterior
|
// Compute the posterior
|
||||||
double logPD = lpd.compute(steadyState, estParams2, deepParams, data, Q, H, presample, info);
|
double logPD = lpd.compute(steadyState, estParams, deepParams, data, Q, H, presample, info);
|
||||||
|
|
||||||
// Cleanups
|
// Cleanups
|
||||||
for (std::vector<EstimatedParameter>::iterator it = estParamsInfo.begin();
|
for (std::vector<EstimatedParameter>::iterator it = estParamsInfo.begin();
|
||||||
|
@ -211,31 +193,80 @@ void
|
||||||
mexFunction(int nlhs, mxArray *plhs[],
|
mexFunction(int nlhs, mxArray *plhs[],
|
||||||
int nrhs, const mxArray *prhs[])
|
int nrhs, const mxArray *prhs[])
|
||||||
{
|
{
|
||||||
if (nrhs != 2 || nlhs != 2)
|
if (nrhs != 7 )
|
||||||
DYN_MEX_FUNC_ERR_MSG_TXT("logposterior: exactly two input arguments and two output arguments are required.");
|
DYN_MEX_FUNC_ERR_MSG_TXT("logposterior: exactly 7 input arguments are required.");
|
||||||
|
|
||||||
// Check and retrieve the arguments
|
if (nlhs > 9 )
|
||||||
|
DYN_MEX_FUNC_ERR_MSG_TXT("logposterior returns 8 output arguments at the most.");
|
||||||
|
|
||||||
|
// Check and retrieve the RHS arguments
|
||||||
|
|
||||||
if (!mxIsDouble(prhs[0]) || mxGetN(prhs[0]) != 1)
|
if (!mxIsDouble(prhs[0]) || mxGetN(prhs[0]) != 1)
|
||||||
DYN_MEX_FUNC_ERR_MSG_TXT("logposterior: First argument must be a column vector of double-precision numbers");
|
DYN_MEX_FUNC_ERR_MSG_TXT("logposterior: First argument must be a column vector of double-precision numbers");
|
||||||
|
|
||||||
VectorConstView estParams(mxGetPr(prhs[0]), mxGetM(prhs[0]), 1);
|
VectorConstView estParams(mxGetPr(prhs[0]), mxGetM(prhs[0]), 1);
|
||||||
|
|
||||||
if (!mxIsDouble(prhs[1]))
|
for (int i = 1; i < 7; ++i)
|
||||||
DYN_MEX_FUNC_ERR_MSG_TXT("logposterior: Second argument must be a matrix of double-precision numbers");
|
if (!mxIsStruct(prhs[i]))
|
||||||
|
{
|
||||||
|
std::stringstream msg;
|
||||||
|
msg << "logposterior: argument " << i+1 << " must be a Matlab structure";
|
||||||
|
DYN_MEX_FUNC_ERR_MSG_TXT(msg.str().c_str());
|
||||||
|
}
|
||||||
|
|
||||||
|
const mxArray *dataset = prhs[1];
|
||||||
|
const mxArray *options_ = prhs[2];
|
||||||
|
const mxArray *M_ = prhs[3];
|
||||||
|
const mxArray *estim_params_ = prhs[4];
|
||||||
|
const mxArray *bayestopt_ = prhs[5];
|
||||||
|
const mxArray *oo_ = prhs[6];
|
||||||
|
|
||||||
|
mxArray *dataset_data = mxGetField(dataset,0,"data");
|
||||||
|
MatrixConstView data(mxGetPr(dataset_data), mxGetM(dataset_data), mxGetN(dataset_data), mxGetM(dataset_data));
|
||||||
|
|
||||||
|
// Creaete LHS arguments
|
||||||
|
|
||||||
|
size_t endo_nbr = (size_t) *mxGetPr(mxGetField(M_, 0, "endo_nbr"));
|
||||||
|
size_t exo_nbr = (size_t) *mxGetPr(mxGetField(M_, 0, "exo_nbr"));
|
||||||
|
size_t param_nbr = (size_t) *mxGetPr(mxGetField(M_, 0, "param_nbr"));
|
||||||
|
size_t varobs_nbr = mxGetM(mxGetField(options_, 0, "varobs"));
|
||||||
|
plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
|
||||||
|
plhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL);
|
||||||
|
plhs[2] = mxCreateDoubleMatrix(endo_nbr, 1, mxREAL);
|
||||||
|
plhs[3] = mxCreateDoubleMatrix(varobs_nbr, 1, mxREAL);
|
||||||
|
plhs[4] = mxCreateDoubleMatrix(1, 1, mxREAL);
|
||||||
|
plhs[5] = mxCreateDoubleMatrix(param_nbr, 1, mxREAL);
|
||||||
|
plhs[6] = mxCreateDoubleMatrix(varobs_nbr, varobs_nbr, mxREAL);
|
||||||
|
plhs[7] = mxCreateDoubleMatrix(exo_nbr, exo_nbr, mxREAL);
|
||||||
|
double *lik = mxGetPr(plhs[0]);
|
||||||
|
double *exit_flag = mxGetPr(plhs[1]);
|
||||||
|
|
||||||
|
VectorView steadyState(mxGetPr(mxGetField(oo_,0,"steady_state")),endo_nbr, 1);
|
||||||
|
VectorView deepParams(mxGetPr(mxGetField(M_, 0, "params")),param_nbr,1);
|
||||||
|
|
||||||
|
MatrixView Q(mxGetPr(mxGetField(M_, 0, "Sigma_e")), exo_nbr, exo_nbr, exo_nbr);
|
||||||
|
|
||||||
|
Matrix H(varobs_nbr,varobs_nbr);
|
||||||
|
const mxArray *H_mx = mxGetField(M_, 0, "H");
|
||||||
|
if (mxGetM(H_mx) == 1 && mxGetN(H_mx) == 1 && *mxGetPr(H_mx) == 0)
|
||||||
|
H.setAll(0.0);
|
||||||
|
else
|
||||||
|
H = MatrixConstView(mxGetPr(H_mx), varobs_nbr, varobs_nbr, varobs_nbr);
|
||||||
|
|
||||||
|
double *trend_coeff = mxGetPr(plhs[3]);
|
||||||
|
double *info_mx = mxGetPr(plhs[4]);
|
||||||
|
|
||||||
MatrixConstView data(mxGetPr(prhs[1]), mxGetM(prhs[1]), mxGetN(prhs[1]), mxGetM(prhs[1]));
|
|
||||||
|
|
||||||
// Compute and return the value
|
// Compute and return the value
|
||||||
try
|
try
|
||||||
{
|
{
|
||||||
double lik = logposterior(estParams, data);
|
int info;
|
||||||
plhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL);
|
*lik = logposterior(estParams, data, options_, M_, estim_params_, bayestopt_, oo_,
|
||||||
*mxGetPr(plhs[1]) = lik;
|
steadyState, trend_coeff, info, deepParams, H, Q);
|
||||||
|
*info_mx = info;
|
||||||
|
*exit_flag = info;
|
||||||
}
|
}
|
||||||
catch (LogposteriorMexErrMsgTxtException e)
|
catch (LogposteriorMexErrMsgTxtException e)
|
||||||
{
|
{
|
||||||
DYN_MEX_FUNC_ERR_MSG_TXT(e.getErrMsg());
|
DYN_MEX_FUNC_ERR_MSG_TXT(e.getErrMsg());
|
||||||
}
|
}
|
||||||
plhs[0] = mxCreateDoubleScalar(0);
|
|
||||||
}
|
}
|
||||||
|
|
|
@ -90,10 +90,3 @@ DynamicModelDLL::~DynamicModelDLL()
|
||||||
#endif
|
#endif
|
||||||
}
|
}
|
||||||
|
|
||||||
void
|
|
||||||
DynamicModelDLL::eval(const Vector &y, const Matrix &x, const Vector &modParams, VectorView &ySteady,
|
|
||||||
Vector &residual, Matrix *g1, Matrix *g2, Matrix *g3) throw (TSException)
|
|
||||||
{
|
|
||||||
Dynamic(y.getData(), x.getData(), 1, modParams.getData(), ySteady.getData(), 0, residual.getData(),
|
|
||||||
g1 == NULL ? NULL : g1->getData(), g2 == NULL ? NULL : g2->getData(), g3 == NULL ? NULL : g3->getData());
|
|
||||||
}
|
|
||||||
|
|
|
@ -59,4 +59,11 @@ public:
|
||||||
//! evaluate Dynamic model DLL
|
//! evaluate Dynamic model DLL
|
||||||
void eval(const Vector &y, const Matrix &x, const Vector ¶ms, VectorView &ySteady,
|
void eval(const Vector &y, const Matrix &x, const Vector ¶ms, VectorView &ySteady,
|
||||||
Vector &residual, Matrix *g1, Matrix *g2, Matrix *g3) throw (TSException);
|
Vector &residual, Matrix *g1, Matrix *g2, Matrix *g3) throw (TSException);
|
||||||
|
template<class VEC>
|
||||||
|
void eval(const Vector &y, const Matrix &x, const VectorView &modParams, VEC &ySteady,
|
||||||
|
Vector &residual, Matrix *g1, Matrix *g2, Matrix *g3) throw (TSException)
|
||||||
|
{
|
||||||
|
Dynamic(y.getData(), x.getData(), 1, modParams.getData(), ySteady.getData(), 0, residual.getData(),
|
||||||
|
g1 == NULL ? NULL : g1->getData(), g2 == NULL ? NULL : g2->getData(), g3 == NULL ? NULL : g3->getData());
|
||||||
|
};
|
||||||
};
|
};
|
||||||
|
|
Loading…
Reference in New Issue