Update, still subject to further testing
git-svn-id: https://www.dynare.org/svn/dynare/trunk@3160 ac1d8469-bf42-47a9-8791-bf33cf982152time-shift
george
parent
98c1291521
commit
d28442f206
|
|
@ -36,7 +36,8 @@ DsgeLikelihood::DsgeLikelihood( Vector& inA_init, GeneralMatrix& inQ, GeneralMa
|
|||
const Vector& INub, const Vector& INlb, const vector<int>&INpshape,
|
||||
const Vector&INp6, const Vector&INp7, const Vector&INp3, const Vector&INp4,
|
||||
Vector& INSteadyState, Vector& INconstant, GeneralParams& INdynareParams, //GeneralParams& parameterDescription,
|
||||
GeneralParams& INdr, GeneralMatrix& INkstate, GeneralMatrix& INghx, GeneralMatrix& INghu
|
||||
GeneralParams& INdr, GeneralMatrix& INkstate, GeneralMatrix& INghx, GeneralMatrix& INghu,
|
||||
GeneralMatrix& INaux, vector<int>&INiv, vector<int>&INic
|
||||
,const int jcols, const char *dfExt)//, KordpDynare& kOrdModel, Approximation& INapprox )
|
||||
:a_init(inA_init), Q(inQ), R(inR), T(inT), Z(inZ), Pstar(inPstar), Pinf(inPinf), H(inH), data(inData), Y(inY),
|
||||
numPeriods(INnumPeriods), numVarobs(inData.numRows()), numTimeObs(inData.numCols()),
|
||||
|
|
@ -47,7 +48,8 @@ DsgeLikelihood::DsgeLikelihood( Vector& inA_init, GeneralMatrix& inQ, GeneralMa
|
|||
param_ub(INub), param_lb(INlb),pshape(INpshape), //pshape(dynareParams.getIntVectorField(string("pshape),
|
||||
p6(INp6), p7(INp7), p3(INp3), p4(INp4), SteadyState(INSteadyState),
|
||||
constant(INconstant), dynareParams(INdynareParams),
|
||||
dr(INdr), kstate(INkstate), ghx(INghx),ghu(INghu)
|
||||
dr(INdr), kstate(INkstate), ghx(INghx),ghu(INghu),
|
||||
aux(INaux), iv(INiv), ic(INic)
|
||||
//, model(kOrdModel), approx(INapprox )
|
||||
{
|
||||
|
||||
|
|
@ -60,24 +62,37 @@ DsgeLikelihood::DsgeLikelihood( Vector& inA_init, GeneralMatrix& inQ, GeneralMa
|
|||
*********/
|
||||
// setting some frequently used common variables that do not need updating
|
||||
//std::vector<double>* vll=new std::vector<double> (nper);
|
||||
vll=new std::vector<double> (numTimeObs);// vector of likelihoods
|
||||
// vll=new std::vector<double> (numTimeObs);// vector of likelihoods
|
||||
vll=new vector<double>(numTimeObs);// vector of likelihoods
|
||||
kalman_algo=(int)dynareParams.getDoubleField(string("kalman_algo"));
|
||||
presampleStart=1+(int)dynareParams.getDoubleField(string("presample"));
|
||||
mode_compute=(int)dr.getDoubleField(string("mode_compute"));
|
||||
mode_compute=(int)dynareParams.getDoubleField(string("mode_compute"));
|
||||
#ifdef DEBUG
|
||||
mexPrintf("mode_compute=%d presampleStart=%d \n", mode_compute,presampleStart);
|
||||
#endif
|
||||
|
||||
// Pepare data for Constructing k-order-perturbation classes
|
||||
//const char *
|
||||
string fname=dynareParams.getStringField(string("fname"));
|
||||
string& fname=dynareParams.getStringField(string("fname"));
|
||||
fName = (char *)fname.c_str();
|
||||
double qz_criterium = dynareParams.getDoubleField(string("qz_criterium"));//qz_criterium = 1+1e-6;
|
||||
int nMax_lag =(int)dynareParams.getDoubleField(string("maximum_lag"));
|
||||
const int nBoth=(int)dr.getDoubleField(string("nboth"));
|
||||
const int nPred = npred-nBoth; // correct nPred for nBoth.
|
||||
//vector<int> *var_order_vp = &order_var;
|
||||
vCov = new TwoDMatrix(Q);
|
||||
#ifdef DEBUG
|
||||
mexPrintf("fName %s qz_criterium=%f nMax_lag=%d nPred=%d :Construction of vCov\n", fName,qz_criterium,nMax_lag,nPred);
|
||||
#endif
|
||||
vCov = new TwoDMatrix(Q);
|
||||
// the lag, current and lead blocks of the jacobian respectively
|
||||
llincidence = new TwoDMatrix (dynareParams.getMatrixField(string("lead_lag_incidence")));
|
||||
#ifdef DEBUG
|
||||
mexPrintf("Construction of casOrdEndoNames\n");
|
||||
#endif
|
||||
charArraySt * casOrdEndoNames=dynareParams.getCharArrayField(string("var_order_endo_names"));
|
||||
#ifdef DEBUG
|
||||
mexPrintf("Construction of endoNamesMX\n");
|
||||
#endif
|
||||
const char **endoNamesMX=(const char ** )casOrdEndoNames->charArrayPtr;
|
||||
|
||||
#ifdef DEBUG
|
||||
|
|
@ -113,7 +128,7 @@ DsgeLikelihood::DsgeLikelihood( Vector& inA_init, GeneralMatrix& inQ, GeneralMa
|
|||
|
||||
// intiate tensor library
|
||||
#ifdef DEBUG
|
||||
mexPrintf("k_order_perturbation: Call tls init\n");
|
||||
mexPrintf("k_order_perturbation: Call tls init order:%d, size: %d\n", order, nstatic+2*nPred+3*nBoth+2*nfwrd+exo_nbr);
|
||||
#endif
|
||||
tls.init(order, nstatic+2*nPred+3*nBoth+2*nfwrd+exo_nbr);
|
||||
|
||||
|
|
@ -203,31 +218,47 @@ DsgeLikelihood::DsgeLikelihood( Vector& inA_init, GeneralMatrix& inQ, GeneralMa
|
|||
|
||||
DsgeLikelihood::~DsgeLikelihood()
|
||||
{
|
||||
|
||||
delete approx;
|
||||
delete model;
|
||||
delete dynamicDLLp;
|
||||
delete journal;
|
||||
delete llincidence;
|
||||
delete vCov;
|
||||
// delete llincidence;
|
||||
// delete vCov;
|
||||
delete vll;
|
||||
/********
|
||||
delete &H;
|
||||
delete &Q;
|
||||
delete &SteadyState;
|
||||
delete &kstate;
|
||||
delete &pshape;
|
||||
******** delete Vectors *********
|
||||
#ifdef DEBUG
|
||||
mexPrintf("delete SS.\n");
|
||||
#endif
|
||||
delete &SteadyState;
|
||||
delete ¶m_ub;
|
||||
delete ¶m_lb;
|
||||
delete &pshape;
|
||||
delete &p6;
|
||||
delete &p7;
|
||||
delete &p3;
|
||||
delete &p4;
|
||||
#ifdef DEBUG
|
||||
mexPrintf("delete params Vectors.\n");
|
||||
#endif
|
||||
delete &xparam1;
|
||||
delete &deepParams;
|
||||
#ifdef DEBUG
|
||||
mexPrintf("delete ghx.\n");
|
||||
#endif
|
||||
**********
|
||||
delete &ghx;
|
||||
delete &ghu;
|
||||
delete &dynareParams;
|
||||
delete &dr;
|
||||
|
||||
delete &aux;
|
||||
delete &iv;
|
||||
delete ⁣
|
||||
**********/
|
||||
}
|
||||
|
||||
double
|
||||
|
|
@ -235,6 +266,7 @@ DsgeLikelihood::CalcLikelihood(Vector& xparams)
|
|||
// runs all routines needed to calculate likelihood
|
||||
{
|
||||
likelihood=0.0;
|
||||
info=0;
|
||||
xparam1=xparams;
|
||||
/*******************************
|
||||
* loop for each sub-sample period
|
||||
|
|
@ -280,18 +312,18 @@ DsgeLikelihood::CalcLikelihood(Vector& xparams)
|
|||
return likelihood;
|
||||
}
|
||||
} // mode compute
|
||||
|
||||
#ifdef DEBUG
|
||||
mexPrintf("Calling of updataeHQparams\n");
|
||||
#endif
|
||||
if(info=updateQHparams()) // updates Q and H matrices and deep parameters
|
||||
{
|
||||
#ifdef DEBUG
|
||||
mexPrintf("Failing of updataeHQparams info =%d\n", info);
|
||||
#endif
|
||||
return likelihood;
|
||||
|
||||
/*****************************************************************************--
|
||||
}
|
||||
/*****************************************************************************--
|
||||
% 2. call model setup & reduction program and pre-filter data
|
||||
******************************************************************************-*/
|
||||
GeneralMatrix& aux = dynareParams.getMatrixField(string("bayestopt_.restrict_aux"));
|
||||
vector<int>&iv= dynareParams.getIntVectorField(string("restrict_var_list"));
|
||||
vector<int>&ic= dynareParams.getIntVectorField(string("bayestopt_.restrict_columns"));
|
||||
|
||||
/*************************************************************
|
||||
// dynare_resolve(() // ... comes here doing:
|
||||
// resol:
|
||||
// check if ys is steady state and calculate one i not
|
||||
|
|
@ -299,36 +331,52 @@ DsgeLikelihood::CalcLikelihood(Vector& xparams)
|
|||
// kalman_transition_matrix(out: A,B, in dr)
|
||||
// IN: bayestopt_.restrict_var_list, bayestopt_.restrict_columns, bayestopt_.restrict_aux, )
|
||||
***************************************************************/
|
||||
#ifdef DEBUG
|
||||
mexPrintf(" *********** Calling dynareResolveDR *********** \n");
|
||||
#endif
|
||||
if (info = dynareResolveDR (iv,ic, aux)) //OUT: [T,R,SteadyState],
|
||||
{
|
||||
#ifdef DEBUG
|
||||
mexPrintf("Failing of dynareResolveDR info =%d\n", info);
|
||||
#endif
|
||||
return likelihood=penalty;
|
||||
}
|
||||
|
||||
/*****************************************************************************--
|
||||
% 2.b pre-filter and detrend data
|
||||
******************************************************************************-*/
|
||||
#ifdef DEBUG
|
||||
mexPrintf("*********** pre-filter and detrend data *********** \n");
|
||||
#endif
|
||||
|
||||
//if options_.noconstant
|
||||
if (dynareParams.getCharField(string("noconstant")))
|
||||
if ((int)dynareParams.getDoubleField(string("noconstant")))
|
||||
constant.zeros();
|
||||
else
|
||||
{
|
||||
//if options_.loglinear
|
||||
if (dynareParams.getCharField(string("loglinear")))
|
||||
if ((int)dynareParams.getDoubleField(string("loglinear")))
|
||||
{
|
||||
for (i =0;i<numVarobs;++i)
|
||||
constant[i] = log(SteadyState[mfys[i]]);
|
||||
constant[i] = log(SteadyState[mfys[i]-1]);
|
||||
}
|
||||
else
|
||||
{
|
||||
for (i =0;i<numVarobs;++i)
|
||||
constant[i] = SteadyState[mfys[i]];
|
||||
constant[i] = SteadyState[mfys[i]-1];
|
||||
}
|
||||
}
|
||||
Vector trend_coeff(numVarobs);
|
||||
//trend = repmat(constant,1,gend);
|
||||
GeneralMatrix constMx(constant.base(),numVarobs,1);
|
||||
#ifdef DEBUG
|
||||
mexPrintf("Calling constMx.repmat numTimeObs=%d\n",numTimeObs);
|
||||
#endif
|
||||
GeneralMatrix&trend = constMx.repmat(1,numTimeObs);
|
||||
|
||||
//if bayestopt_.with_trend
|
||||
if (dynareParams.getCharField(string("with_trend")))
|
||||
/************************************
|
||||
if ((int)dynareParams.getDoubleField(string("with_trend")))
|
||||
{
|
||||
trend_coeff.zeros();
|
||||
// GeneralMatrix& mt = dynareParams.getMatrixField(string("trend_coeffs"));
|
||||
|
|
@ -347,15 +395,22 @@ DsgeLikelihood::CalcLikelihood(Vector& xparams)
|
|||
|
||||
trend.add(1,trend_coefMx);
|
||||
}
|
||||
*************************************/
|
||||
presampleStart =(int) dynareParams.getDoubleField(string("presample"))+1;
|
||||
int no_missing_data_flag = (number_of_observations==numTimeObs*numVarobs);
|
||||
//Y =data-trend;
|
||||
Y=data;
|
||||
#ifdef DEBUG
|
||||
mexPrintf("Calling Y.add( trend) in GeneralMatrices\n");
|
||||
#endif
|
||||
Y.add(-1,trend);
|
||||
|
||||
/*****************************************************************************
|
||||
% 3. Initial condition of the Kalman filter
|
||||
*******************************************************************************/
|
||||
#ifdef DEBUG
|
||||
mexPrintf("Calling InitiateKalmanMatrices\n");
|
||||
#endif
|
||||
if( InitiateKalmanMatrices())
|
||||
return likelihood=penalty;
|
||||
|
||||
|
|
@ -364,14 +419,16 @@ DsgeLikelihood::CalcLikelihood(Vector& xparams)
|
|||
// choose and run KF to get likelihood fval
|
||||
*****************************************************************************/
|
||||
likelihood+=KalmanFilter(0.000001, false);// calls Kalman
|
||||
|
||||
|
||||
/****************************************************************************
|
||||
// Adds prior if necessary
|
||||
****************************************************************************/
|
||||
//likelihood-= priordens(xparam1,pshape,p6,p7,p3,p4);//fval = (likelihood-lnprior);
|
||||
//options_.kalman_algo = kalman_algo;
|
||||
} // end sub-sample loop
|
||||
#ifdef DEBUG
|
||||
mexPrintf("End of CallcLlikelihood returning likelihood=%f\n", likelihood);
|
||||
#endif
|
||||
return likelihood;
|
||||
}
|
||||
|
||||
/**************************************************
|
||||
|
|
@ -386,11 +443,23 @@ int
|
|||
DsgeLikelihood::updateQHparams()// updates Q and H matrices and deep parameters
|
||||
{
|
||||
int i=0, offset=0, nv=0, k, k1, k2, info=0;
|
||||
#ifdef DEBUG
|
||||
mexPrintf("Setting of Q \n");
|
||||
#endif
|
||||
delete &Q;
|
||||
Q = dynareParams.getMatrixField(string("Sigma_e"));
|
||||
#ifdef DEBUG
|
||||
mexPrintf("Setting of H \n");
|
||||
#endif
|
||||
delete &H;
|
||||
H = dynareParams.getMatrixField(string("H"));
|
||||
nv=(int)dynareParams.getDoubleField(string("nvx"));
|
||||
if(nv)
|
||||
{
|
||||
#ifdef DEBUG
|
||||
mexPrintf("Setting of Q var_exo\n");
|
||||
#endif
|
||||
// if(&estvx) delete &estvx;
|
||||
GeneralMatrix&estvx=dynareParams.getMatrixField(string("var_exo"));
|
||||
for (i=0;i<nv;++i)
|
||||
{
|
||||
|
|
@ -398,11 +467,14 @@ DsgeLikelihood::updateQHparams()// updates Q and H matrices and deep parameters
|
|||
Q.get(k,k) = xparam1[i]*xparam1[i];
|
||||
}
|
||||
offset = nv;
|
||||
delete &estvx;
|
||||
}
|
||||
|
||||
nv=(int)dynareParams.getDoubleField(string("nvn"));
|
||||
if(nv)
|
||||
{
|
||||
#ifdef DEBUG
|
||||
mexPrintf("Setting of H var_endo\n");
|
||||
#endif
|
||||
GeneralMatrix&estvn=dynareParams.getMatrixField(string("var_endo"));
|
||||
for (i=0;i<nv;++i)
|
||||
{
|
||||
|
|
@ -410,13 +482,16 @@ DsgeLikelihood::updateQHparams()// updates Q and H matrices and deep parameters
|
|||
H.get(k,k) = xparam1[i+offset]*xparam1[i+offset];
|
||||
}
|
||||
offset += nv;
|
||||
delete &estvn;
|
||||
}
|
||||
|
||||
//if estim_params_.ncx
|
||||
//for i=1:estim_params_.ncx
|
||||
nv=(int)dynareParams.getDoubleField(string("ncx"));
|
||||
if(nv)
|
||||
{
|
||||
#ifdef DEBUG
|
||||
mexPrintf("Setting of Q corrx\n");
|
||||
#endif
|
||||
GeneralMatrix&corrx=dynareParams.getMatrixField(string("corrx"));
|
||||
for (i=0;i<nv;++i)
|
||||
{
|
||||
|
|
@ -426,6 +501,7 @@ DsgeLikelihood::updateQHparams()// updates Q and H matrices and deep parameters
|
|||
Q.get(k2,k1) = Q.get(k1,k2);
|
||||
}
|
||||
// [CholQ,testQ] = chol(Q);
|
||||
delete &corrx;
|
||||
int testQ=0;
|
||||
try
|
||||
{
|
||||
|
|
@ -472,6 +548,9 @@ DsgeLikelihood::updateQHparams()// updates Q and H matrices and deep parameters
|
|||
nv=(int)dynareParams.getDoubleField(string("ncn"));
|
||||
if(nv)
|
||||
{
|
||||
#ifdef DEBUG
|
||||
mexPrintf("Setting of H corrn\n");
|
||||
#endif
|
||||
GeneralMatrix&corrn=dynareParams.getMatrixField(string("corrn"));
|
||||
vector<int>&lgyidx2varobs= dynareParams.getIntVectorField(string("lgyidx2varobs"));
|
||||
for (i=0;i<nv;++i)
|
||||
|
|
@ -485,7 +564,8 @@ DsgeLikelihood::updateQHparams()// updates Q and H matrices and deep parameters
|
|||
H.get(k1,k2) = xparam1[i+offset]*sqrt(H.get(k1,k1)*H.get(k2,k2));
|
||||
H.get(k2,k1) = H.get(k1,k2);
|
||||
}
|
||||
|
||||
delete &corrn;
|
||||
|
||||
//[CholH,testH] = chol(H);
|
||||
int testH=0;
|
||||
try
|
||||
|
|
@ -535,6 +615,9 @@ DsgeLikelihood::updateQHparams()// updates Q and H matrices and deep parameters
|
|||
else
|
||||
TS_RAISE("Inssuficient length of the xparam1 parameters vector");
|
||||
}
|
||||
#ifdef DEBUG
|
||||
mexPrintf("End of Setting of HQ params\n");
|
||||
#endif
|
||||
|
||||
/**********
|
||||
M_.Sigma_e = Q;
|
||||
|
|
@ -551,37 +634,67 @@ DsgeLikelihood::InitiateKalmanMatrices()
|
|||
{
|
||||
int np = T.numRows();// size(T,1);
|
||||
double lyapunov_tol=dynareParams.getDoubleField(string("lyapunov_tol"));
|
||||
int lik_init=(int)dynareParams.getDoubleField(string("lik_init"));
|
||||
// int lik_init=(int)dynareParams.getDoubleField(string("lik_init"));
|
||||
//if options_.lik_init == 1 % Kalman filter
|
||||
GeneralMatrix RQRt(R,Q); // R*Q
|
||||
RQRt.multRightTrans(R); // R*Q*Rt
|
||||
// GeneralMatrix RQRt(R,Q); // R*Q
|
||||
GeneralMatrix RQ(R,Q); // R*Q
|
||||
#ifdef DEBUG
|
||||
mexPrintf("Calling RQRt.multRightTrans(R)\n");
|
||||
#endif
|
||||
// GeneralMatrix::md_length=RQRt.numRows();
|
||||
//RQRt.md_length=RQRt.numRows();
|
||||
|
||||
// RQRt.multRightTrans(R); // R*Q*Rt
|
||||
// GeneralMatrix RQRt(np,np);
|
||||
// RQRt.zeros();
|
||||
// RQRt.multAndAdd(RQ,R, "T", 1.0);
|
||||
GeneralMatrix RQRt(RQ,R, "T");
|
||||
#ifdef DEBUG
|
||||
for (int i=0;i<RQRt.numRows();++i)
|
||||
{
|
||||
for (int j=0;j<RQRt.numCols();++j)
|
||||
mexPrintf(" %f ", RQRt.get(i,j));
|
||||
mexPrintf("\n");
|
||||
}
|
||||
#endif
|
||||
GeneralMatrix Ptmp(np,np);
|
||||
//Pstar = lyapunov_symm (T,R*Q*R',options_.qz_criterium,options_.lyapunov_complex_threshold)
|
||||
#ifdef DEBUG
|
||||
mexPrintf("Calling disclyap_fast to initialise Pstar:\n");
|
||||
#endif
|
||||
disclyap_fast(T,RQRt,Ptmp, lyapunov_tol, 1); // 1 to check chol
|
||||
Pstar=Ptmp;
|
||||
//Pinf=[]
|
||||
//Pinf = *(new GeneralMatrix(np,np));
|
||||
#ifdef DEBUG
|
||||
Pstar.print();
|
||||
mexPrintf("Initialise Pinf\n");
|
||||
#endif
|
||||
Pinf.zeros();
|
||||
|
||||
//Z= zeros(size(data,1), size(T,2))
|
||||
GeneralMatrix Ztmp=*(new GeneralMatrix(numVarobs,np));
|
||||
Ztmp.zeros();
|
||||
|
||||
//a=zeros(size(T,1),1);
|
||||
Vector atmp(np);
|
||||
atmp.zeros();
|
||||
a_init=atmp;
|
||||
|
||||
//Vector atmp(np);
|
||||
//atmp.zeros();
|
||||
//a_init=atmp;
|
||||
a_init.zeros();
|
||||
//if (lik_init == 2)// Old Diffuse Kalman filter
|
||||
// Pstar = options_.Harvey_scale_factor*eye(np);
|
||||
//Pinf = [];
|
||||
//else if (lik_init == 3) // Diffuse Kalman filter
|
||||
// else ...
|
||||
|
||||
for (int i = 0;i<numVarobs;++i)
|
||||
Ztmp.get(i,mf[i])=1.0;
|
||||
Z=Ztmp;
|
||||
delete &Ztmp;
|
||||
#ifdef DEBUG
|
||||
mexPrintf("Initialise Z\n");
|
||||
#endif
|
||||
//Z= zeros(size(data,1), size(T,2))
|
||||
//Z.zeros();
|
||||
//GeneralMatrix Ztmp=*(new GeneralMatrix(numVarobs,np));
|
||||
//for (int i = 0;i<numVarobs;++i)
|
||||
// Ztmp.get(i,mf[i]-1)=1.0;
|
||||
//Z=Ztmp;
|
||||
//delete &Ztmp;
|
||||
//for (int i = 0;i<numVarobs;++i)
|
||||
// Z.get(i,mf[i]-1)=1.0;
|
||||
}
|
||||
|
||||
|
||||
|
|
@ -593,15 +706,16 @@ double
|
|||
DsgeLikelihood::KalmanFilter(double riccatiTol=0.000001,bool uni = false)
|
||||
{
|
||||
bool diffuse=false;
|
||||
double loglik;
|
||||
|
||||
try
|
||||
{
|
||||
// make input matrices
|
||||
int start = presampleStart;
|
||||
|
||||
int nper=Y.numCols();
|
||||
#ifdef DEBUG
|
||||
mexPrintf("kalman_filter: periods=%d start=%d, a.length=%d, uni=%d diffuse=%d\n", nper, start,a_init.length(), uni, diffuse);
|
||||
Pstar.print(); Z.print(); H.print(); T.print(); R.print(); Q.print(); Y.print();
|
||||
#endif
|
||||
|
||||
// make storage for output
|
||||
|
|
@ -625,20 +739,21 @@ DsgeLikelihood::KalmanFilter(double riccatiTol=0.000001,bool uni = false)
|
|||
if (uni)
|
||||
{
|
||||
KalmanUniTask kut(kt);
|
||||
likelihood = kut.filter(per, d, (start-1), vll);
|
||||
loglik = kut.filter(per, d, (start-1), vll);
|
||||
per = per / Y.numRows();
|
||||
d = d / Y.numRows();
|
||||
}
|
||||
else
|
||||
{
|
||||
#ifdef TIMING_LOOP
|
||||
mexPrintf("kalman_filter: starting 1000 loops\n");
|
||||
for (int tt=0;tt<1000;++tt)
|
||||
{
|
||||
#endif
|
||||
likelihood = kt.filter(per, d, (start-1), vll);
|
||||
loglik = kt.filter(per, d, (start-1), vll);
|
||||
#ifdef TIMING_LOOP
|
||||
}
|
||||
mexPrintf("kalman_filter: finished 1000 loops");
|
||||
mexPrintf("kalman_filter: finished 1000 loops\n");
|
||||
#endif
|
||||
}
|
||||
}
|
||||
|
|
@ -647,16 +762,17 @@ DsgeLikelihood::KalmanFilter(double riccatiTol=0.000001,bool uni = false)
|
|||
init = new StateInit(Pstar, a_init);
|
||||
BasicKalmanTask bkt(Y, Z, H, T, R, Q, *init, riccatiTol);
|
||||
#ifdef TIMING_LOOP
|
||||
mexPrintf("kalman_filter: starting 1000 loops\n");
|
||||
for (int tt=0;tt<1000;++tt)
|
||||
{
|
||||
#endif
|
||||
likelihood = bkt.filter( per, d, (start-1), vll);
|
||||
loglik = bkt.filter( per, d, (start-1), vll);
|
||||
#ifdef DEBUG
|
||||
mexPrintf("Basickalman_filter: likelihood=%f \n", likelihood);
|
||||
mexPrintf("Basickalman_filter: loglik=%f \n", loglik);
|
||||
#endif
|
||||
#ifdef TIMING_LOOP
|
||||
}
|
||||
mexPrintf("Basickalman_filter: finished 1000 loops");
|
||||
mexPrintf("Basickalman_filter: finished 1000 loops\n");
|
||||
#endif
|
||||
|
||||
}
|
||||
|
|
@ -673,6 +789,6 @@ DsgeLikelihood::KalmanFilter(double riccatiTol=0.000001,bool uni = false)
|
|||
e.printMessage(mes, 299);
|
||||
mexErrMsgTxt(mes);
|
||||
}
|
||||
return likelihood;
|
||||
return loglik;
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -80,7 +80,11 @@ class DsgeLikelihood
|
|||
GeneralParams& dynareParams;
|
||||
//GeneralParams& parameterDescription;
|
||||
GeneralParams& dr;
|
||||
|
||||
|
||||
GeneralMatrix& aux; //= dynareParams.getMatrixField(string("restrict_aux"));
|
||||
vector<int>&iv; //= dynareParams.getIntVectorField(string("restrict_var_list"));
|
||||
vector<int>⁣ //= dynareParams.getIntVectorField(string("restrict_columns"));
|
||||
|
||||
GeneralMatrix& kstate;
|
||||
GeneralMatrix& ghx;
|
||||
GeneralMatrix& ghu;
|
||||
|
|
@ -97,7 +101,7 @@ class DsgeLikelihood
|
|||
//KalmanUniTask ukt;// univariate
|
||||
// member functions
|
||||
MexStructParam& SetDRModel(MexStructParam¶ms);
|
||||
void disclyap_fast(const GeneralMatrix &G, const GeneralMatrix & V, GeneralMatrix &X, double tol = 1e-16, int flag_ch=0);
|
||||
// void disclyap_fast(const GeneralMatrix &G, const GeneralMatrix & V, GeneralMatrix &X, double tol = 1e-16, int flag_ch=0);
|
||||
GeneralMatrix& SolveDRkOrderPert();//calls k-order pert or whatever;
|
||||
int dynareResolveDR(vector<int>&iv,vector<int>&ic,GeneralMatrix& aux); // returns int info, ys, and TT and RR Decision Rule
|
||||
int SolveDRModel(const int endo_nbr, const int exo_nbr, const int nstatic, const int npred, int nfwrd);//int dr1(); // returns int info and updated dr
|
||||
|
|
@ -119,7 +123,8 @@ class DsgeLikelihood
|
|||
const Vector&INp6, const Vector&INp7, const Vector&INp3, const Vector&INp4,
|
||||
Vector& INSteadyState, Vector& INconstant, GeneralParams& INdynareParams,
|
||||
//GeneralParams& parameterDescription,
|
||||
GeneralParams& INdr, GeneralMatrix& INkstate, GeneralMatrix& INghx, GeneralMatrix& INghu
|
||||
GeneralParams& INdr, GeneralMatrix& INkstate, GeneralMatrix& INghx, GeneralMatrix& INghu,
|
||||
GeneralMatrix& aux, vector<int>&iv, vector<int>&ic
|
||||
,const int jcols, const char *dfExt); //, KordpDynare& inModel, Approximation& INapprox );
|
||||
DsgeLikelihood( const Vector¶ms,const GeneralMatrix&data, const vector<int>& data_index, const int gend,
|
||||
const int number_of_observations, const bool no_more_missing_observations);//, KordpDynare& model ); // constructor, and
|
||||
|
|
|
|||
|
|
@ -24,24 +24,38 @@ MexStruct::MexStruct( const int numparstruct):
|
|||
numParStruct(numparstruct), structName(string(""))
|
||||
{
|
||||
// get Dynare mexSturcture pointers and store them locally
|
||||
#ifdef DEBUG
|
||||
mexPrintf("MexStruct reserve=%d \n", numParStruct);
|
||||
#endif
|
||||
parStruct.reserve(numParStruct);
|
||||
parStructBase.reserve(numParStructs);
|
||||
for (int i=0;i<numParStructs;++i)
|
||||
parStructBase.reserve(numParStruct);
|
||||
for (int i=0;i<numParStruct;++i)
|
||||
{
|
||||
parStruct[i]=mexGetVariablePtr("caller", DynareParamStructsNm[i]);
|
||||
parStructBase[i]=caller;
|
||||
#ifdef DEBUG
|
||||
mexPrintf("MexStruct to insert i=%d parStructNm[i]=%s using base[i]=%d %s\n", i, DynareParamStructsNm[i], parStructBase[i], mexBase[parStructBase[i]]);
|
||||
#endif
|
||||
// get field names into the map:
|
||||
pair <map<string, int>::iterator, bool> ret;
|
||||
int j=0;
|
||||
const char* field;
|
||||
while ((field=mxGetFieldNameByNumber(parStruct[i],j))!=NULL)
|
||||
{
|
||||
#ifdef DEBUG
|
||||
mexPrintf("MexStruct insert field= %s\n", field);
|
||||
#endif
|
||||
ret=parNamStructMap.insert(make_pair(string(field),i));
|
||||
if (!ret.second)
|
||||
mexErrMsgTxt("MexStruct: Failed to insert param \n");
|
||||
mexPrintf("MexStruct failed to insert field %s from struct %d %s using base %d %s\n"
|
||||
, field, i, DynareParamStructsNm[i], parStructBase[i], mexBase[parStructBase[i]]);
|
||||
// mexErrMsgTxt("MexStruct: Failed to insert param \n");
|
||||
j++;
|
||||
}
|
||||
}
|
||||
#ifdef DEBUG
|
||||
mexPrintf("MexStruct insert finished \n");
|
||||
#endif
|
||||
}
|
||||
|
||||
MexStructParam&
|
||||
|
|
@ -119,7 +133,7 @@ MexStruct::getIntVectorField(const string& field)
|
|||
{
|
||||
(*vars)[v] = (int)(*(dparams++)); //[v];
|
||||
#ifdef DEBUG
|
||||
mexPrintf("vars)[%d] = %d .\n", v, (*vars)[v]);
|
||||
mexPrintf("%s[%d]=%d.\n", field.c_str() , v, (*vars)[v]);
|
||||
#endif
|
||||
}
|
||||
|
||||
|
|
@ -285,7 +299,7 @@ MexStructParam::getIntVectorField(const string& field)
|
|||
{
|
||||
(*vars)[v] = (int)(*(dparams++)); //[v];
|
||||
#ifdef DEBUG
|
||||
mexPrintf("vars)[%d] = %d .\n", v, (*vars)[v]);
|
||||
mexPrintf("%s[%d]=%d.\n", field.c_str(), v, (*vars)[v]);
|
||||
#endif
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -73,7 +73,10 @@ class MexStruct :public virtual GeneralParams
|
|||
{
|
||||
map<string, int>::iterator it=parNamStructMap.find(field);
|
||||
if (it==parNamStructMap.end())
|
||||
throw(SYLV_MES_EXCEPTION("no parameter with such name"));
|
||||
{
|
||||
mexPrintf("getMxField:no parameter with such name");
|
||||
throw(SYLV_MES_EXCEPTION("getMxField:no parameter with such name"));
|
||||
}
|
||||
return mxGetField(parStruct[it->second], 0, field.c_str() );
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -111,6 +111,9 @@ DsgeLikelihood::SolveDRModel(const int endo_nbr, const int exo_nbr, const int ns
|
|||
//dr.ghu = repmat(1./dr.ys(k1),1,size(dr.ghu,2)).*dr.ghu;
|
||||
GeneralMatrix&repInvSSu=mInvOrdSS.repmat(1,ghu.numCols());
|
||||
ghu.multElements(repInvSSu);
|
||||
delete &repk1klagSSt;
|
||||
delete &repInvSSu;
|
||||
delete &repInvSSx;
|
||||
};//end if
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -59,24 +59,38 @@ DsgeLikelihood::dynareResolveDR(vector<int>&iv,vector<int>&ic,GeneralMatrix& aux
|
|||
// check if ys is steady state and calcluate new one if not
|
||||
// testing for steadystate file: To Be Icluded at a later stage
|
||||
|
||||
|
||||
int info = SolveDRModel(endo_nbr,exo_nbr,nstatic, npred, nfwrd); //formerly known as dr1 in Matalb Dynare, i.e.
|
||||
#ifdef DEBUG
|
||||
mexPrintf("Calling SolveDRModel\n");
|
||||
#endif
|
||||
// [dr,info,M_,options_,oo_] = dr1(dr,check_flag,M_,options_,oo_);
|
||||
|
||||
// End of resol:
|
||||
// now rest of dynare_resolve:
|
||||
int DRinfo = SolveDRModel(endo_nbr,exo_nbr,nstatic, npred, nfwrd); //formerly known as dr1 in Matalb Dynare, i.e.
|
||||
if (DRinfo) return DRinfo;
|
||||
|
||||
// End of resol: now rest of dynare_resolve:
|
||||
|
||||
// if nargin == 0
|
||||
// if (iv.size()==0)
|
||||
#ifdef DEBUG
|
||||
mexPrintf(" if iv==NULL\n");
|
||||
#endif
|
||||
if (&iv==NULL)
|
||||
{
|
||||
#ifdef DEBUG
|
||||
mexPrintf(" construct iv\n");
|
||||
#endif
|
||||
//iv = (1:endo_nbr)';
|
||||
for (int i=1;i<=endo_nbr;++i)
|
||||
iv.push_back(i);//= (1:endo_nbr)';
|
||||
}
|
||||
// if (ic.size()==0)
|
||||
#ifdef DEBUG
|
||||
mexPrintf(" if ic==NULL\n");
|
||||
#endif
|
||||
if (&ic==NULL)
|
||||
{
|
||||
#ifdef DEBUG
|
||||
mexPrintf(" construct ic\n");
|
||||
#endif
|
||||
//ic = [ nstatic+(1:npred) endo_nbr+(1:size(oo_.dr.ghx,2)-npred) ]';
|
||||
//ic(npred+ghx.numCols());
|
||||
for(int i=0;i<npred;++i)
|
||||
|
|
@ -87,8 +101,14 @@ DsgeLikelihood::dynareResolveDR(vector<int>&iv,vector<int>&ic,GeneralMatrix& aux
|
|||
ic.push_back(j+endo_nbr+1);
|
||||
}
|
||||
|
||||
#ifdef DEBUG
|
||||
mexPrintf(" if aux==NULL\n");
|
||||
#endif
|
||||
if (&aux==NULL)
|
||||
{
|
||||
#ifdef DEBUG
|
||||
mexPrintf(" construct aux\n");
|
||||
#endif
|
||||
int i;
|
||||
aux =dr.getMatrixField(string("transition_auxiliary_variables"));
|
||||
//k = find(aux(:,2) > npred);
|
||||
|
|
@ -106,52 +126,68 @@ DsgeLikelihood::dynareResolveDR(vector<int>&iv,vector<int>&ic,GeneralMatrix& aux
|
|||
aux.get(k[i]-1,1) +=nfwrd;
|
||||
}//end if
|
||||
|
||||
#ifdef DEBUG
|
||||
mexPrintf("[A,B] = kalman_transition_matrix\n");
|
||||
#endif
|
||||
|
||||
// here is content of [A,B] = kalman_transition_matrix(oo_.dr,iv,ic,aux,M_.exo_nbr);
|
||||
{
|
||||
int n_iv = iv.size();//length(iv);
|
||||
int n_ir1 = aux.numCols();// size(aux,1);
|
||||
int nr = n_iv + n_ir1;
|
||||
// here is the content of [T R]=[A,B] = kalman_transition_matrix(oo_.dr,iv,ic,aux,M_.exo_nbr);
|
||||
|
||||
int n_iv = iv.size();//length(iv);
|
||||
int n_ir1 = aux.numRows();// size(aux,1);
|
||||
int nr = n_iv + n_ir1;
|
||||
|
||||
GeneralMatrix A=*(new GeneralMatrix (nr,nr));
|
||||
A.zeros();
|
||||
GeneralMatrix B=*(new GeneralMatrix(nr,exo_nbr));
|
||||
B.zeros();
|
||||
|
||||
vector<int>i_n_iv(n_iv);
|
||||
for (int i=0;i<n_iv;++i)
|
||||
i_n_iv[i]=i+1;//= (1:n_iv)';
|
||||
|
||||
|
||||
//A(i_n_iv,ic) = dr.ghx(iv,:);
|
||||
A.AssignByVectors (i_n_iv,ic, ghx, iv, nullVec);//dr.ghx(iv,:);
|
||||
// GeneralMatrix A=*(new GeneralMatrix (nr,nr));
|
||||
// A.zeros();
|
||||
// GeneralMatrix B=*(new GeneralMatrix(nr,exo_nbr));
|
||||
// B.zeros();
|
||||
T.zeros();
|
||||
R.zeros();
|
||||
vector<int>i_n_iv(n_iv);
|
||||
for (int i=0;i<n_iv;++i)
|
||||
i_n_iv[i]=i+1;//= (1:n_iv)';
|
||||
|
||||
#ifdef DEBUG
|
||||
mexPrintf("T=A assign by vec\n");
|
||||
#endif
|
||||
//A(i_n_iv,ic) = dr.ghx(iv,:);
|
||||
//A.AssignByVectors (i_n_iv,ic, ghx, iv, nullVec);//dr.ghx(iv,:);
|
||||
T.AssignByVectors (i_n_iv,ic, ghx, iv, nullVec);//dr.ghx(iv,:);
|
||||
#ifdef DEBUG
|
||||
mexPrintf("Completed T=A assign by vec\n");
|
||||
#endif
|
||||
if (n_ir1 > 0)
|
||||
{
|
||||
//A(n_iv+1:end,:) = sparse(aux(:,1),aux(:,2),ones(n_ir1,1),n_ir1,nr);
|
||||
|
||||
if (n_ir1!=aux.numRows())// throw error
|
||||
throw SYLV_MES_EXCEPTION("Wrong dimensions for aux matrix.");
|
||||
|
||||
GeneralMatrix sparse(n_ir1,nr);
|
||||
for (int i=0;i<n_ir1;++i)
|
||||
sparse.get((int)aux.get(i,0)-1,(int)aux.get(i,1)-1)=1;
|
||||
|
||||
/* vector<int>span2end(A.numRows()-n_iv);
|
||||
for (int i=0;i<A.numRows()-n_iv;++i)
|
||||
{
|
||||
//A(n_iv+1:end,:) = sparse(aux(:,1),aux(:,2),ones(n_ir1,1),n_ir1,nr);
|
||||
|
||||
#ifdef DEBUG
|
||||
mexPrintf("Create sparse\n");
|
||||
#endif
|
||||
GeneralMatrix sparse(n_ir1,nr);
|
||||
for (int i=0;i<n_ir1;++i)
|
||||
sparse.get((int)aux.get(i,0)-1,(int)aux.get(i,1)-1)=1;
|
||||
|
||||
/* vector<int>span2end(A.numRows()-n_iv);
|
||||
for (int i=0;i<A.numRows()-n_iv;++i)
|
||||
span2end[i]=i+n_iv+1;
|
||||
*/
|
||||
A.place(sparse,n_iv,0);
|
||||
}
|
||||
T=A;
|
||||
|
||||
//B(i_n_iv,:) = dr.ghu(iv,:);
|
||||
GeneralMatrix& ghu=dr.getMatrixField(string("ghu"));
|
||||
B.AssignByVectors (i_n_iv, nullVec, ghu, iv, nullVec);
|
||||
|
||||
R=B;
|
||||
|
||||
A.place(sparse,n_iv,0); and T=A;
|
||||
*/
|
||||
#ifdef DEBUG
|
||||
mexPrintf("T.place (sparse,n_iv,0)\n");
|
||||
#endif
|
||||
T.place(sparse,n_iv,0);
|
||||
}
|
||||
|
||||
#ifdef DEBUG
|
||||
mexPrintf("R=B assign by vec R rows %d cols %d \n", R.numRows(), R.numCols());
|
||||
#endif
|
||||
//B(i_n_iv,:) = dr.ghu(iv,:); and R=B;
|
||||
R.AssignByVectors(i_n_iv, nullVec, ghu, iv, nullVec);
|
||||
#ifdef DEBUG
|
||||
mexPrintf("Completed R=B assign by vec\n");
|
||||
#endif
|
||||
|
||||
// ys = oo_.dr.ys;
|
||||
GeneralMatrix&ysmx = dr.getMatrixField(string("ys"));
|
||||
SteadyState=ysmx.getData();
|
||||
// GeneralMatrix&ysmx = dr.getMatrixField(string("ys"));
|
||||
// SteadyState=ysmx.getData();
|
||||
return DRinfo;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -24,8 +24,9 @@
|
|||
/******************************************
|
||||
* mexFunction: Matlab Inerface point and the main application driver
|
||||
* for DsgeLikelihood
|
||||
*****************************************************************%function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data,data_index,number_of_observations,no_more_missing_observations)
|
||||
% function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data,data_index,number_of_observations,no_more_missing_observations)
|
||||
*****************************************************************
|
||||
% function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data,data_index,
|
||||
% number_of_observations,no_more_missing_observations)
|
||||
% Evaluates the posterior kernel of a dsge model.
|
||||
%
|
||||
% INPUTS
|
||||
|
|
@ -49,9 +50,10 @@
|
|||
*****************************************************************/
|
||||
#include "DsgeLikelihood.h"
|
||||
|
||||
#ifdef MATLAB
|
||||
#include "mexutils.h"
|
||||
#endif
|
||||
|
||||
extern const char *DynareParamStructsNm []={"M_", "oo_", "options_", "bayestopt_", "estim_params_", "dr"};
|
||||
extern const char* mexBase[]={"base", "caller", "global"};
|
||||
|
||||
extern "C" {
|
||||
|
||||
|
|
@ -83,24 +85,39 @@ mexFunction(int nlhs, mxArray *plhs[],
|
|||
}
|
||||
|
||||
#ifdef DEBUG
|
||||
mexPrintf("k_order_perturbation: mexExt=%s.\n", dfExt);
|
||||
mexPrintf("estimation: mexExt=%s.\n", dfExt);
|
||||
#endif
|
||||
/***********
|
||||
***************/
|
||||
int numPeriods=1;
|
||||
|
||||
//MexStruct dynareParams();
|
||||
GeneralParams& dynareParams=*(new MexStruct());
|
||||
//MexStructParam& dr=dynareParams.getStructField("dr");
|
||||
GeneralParams& dr=dynareParams.getStructField("dr");
|
||||
|
||||
MexStruct& dynareParams=*(new MexStruct(numParStructs));
|
||||
#ifdef DEBUG
|
||||
mexPrintf("getting dr\n");
|
||||
#endif
|
||||
MexStructParam& dr=dynareParams.getMexStructField("dr");
|
||||
vector<int>&mfys=dynareParams.getIntVectorField("mfys");
|
||||
vector<int>&mf=dynareParams.getIntVectorField("mf1");
|
||||
int numVarobs=data.numRows();
|
||||
#ifdef DEBUG
|
||||
mexPrintf("getting SS\n");
|
||||
#endif
|
||||
Vector& SteadyState=dr.getDoubleVectorField(string("ys"));
|
||||
#ifdef DEBUG
|
||||
int gg;
|
||||
for ( gg=0;gg<SteadyState.length();++gg)
|
||||
mexPrintf("SteadyState %d = %f\n", gg, SteadyState[gg]);
|
||||
#endif
|
||||
|
||||
int numVarobs=data.numRows();
|
||||
|
||||
Vector constant(numVarobs);//=*(new Vector(numVarobs));// = *(new Vector(nobs));
|
||||
GeneralMatrix&kstate = dr.getMatrixField(string("kstate"));
|
||||
vector<int>&order_var = dr.getIntVectorField(string("order_var"));
|
||||
#ifdef DEBUG
|
||||
for ( gg=0;gg<order_var.size();++gg)
|
||||
mexPrintf("order_var %d = %d\n", gg, order_var[gg]);
|
||||
#endif
|
||||
int order=(int)dynareParams.getDoubleField(string("order"));
|
||||
int endo_nbr = (int)dynareParams.getDoubleField(string("endo_nbr"));
|
||||
int exo_nbr = (int)dynareParams.getDoubleField(string("exo_nbr"));
|
||||
|
|
@ -109,6 +126,10 @@ mexFunction(int nlhs, mxArray *plhs[],
|
|||
int nfwrd = (int)dr.getDoubleField(string("nfwrd"));
|
||||
Vector& ub=dynareParams.getDoubleVectorField(string("ub"));
|
||||
Vector& lb=dynareParams.getDoubleVectorField(string("lb"));
|
||||
#ifdef DEBUG
|
||||
for ( gg=0;gg<lb.length();++gg)
|
||||
mexPrintf("lb %d = %f\n", gg, lb[gg]);
|
||||
#endif
|
||||
int num_dp=(int)dynareParams.getDoubleField(string("np"));// no of deep params
|
||||
Vector& deepParams=*(new Vector(num_dp));
|
||||
vector<int>&pshape=dynareParams.getIntVectorField(string("pshape"));
|
||||
|
|
@ -121,36 +142,69 @@ mexFunction(int nlhs, mxArray *plhs[],
|
|||
int nsPred=(int)dr.getDoubleField(string("nspred"));
|
||||
int nsForw=(int)dr.getDoubleField(string("nsfwrd"));
|
||||
const int jcols = exo_nbr+endo_nbr+nsPred+nsForw;
|
||||
#ifdef DEBUG
|
||||
mexPrintf("jcols = %d, exo_nbr=%d\n", jcols, exo_nbr);
|
||||
#endif
|
||||
|
||||
GeneralMatrix& aux= dr.getMatrixField(string("transition_auxiliary_variables"));
|
||||
int nr=endo_nbr+aux.numRows();
|
||||
Vector& a_init=*(new Vector(numVarobs));
|
||||
GeneralMatrix& aux = dynareParams.getMatrixField(string("restrict_aux"));
|
||||
vector<int>&iv= dynareParams.getIntVectorField(string("restrict_var_list"));
|
||||
vector<int>&ic= dynareParams.getIntVectorField(string("restrict_columns"));
|
||||
|
||||
int nr=nsPred+aux.numRows();
|
||||
Vector& a_init=*(new Vector(nr));
|
||||
a_init.zeros();
|
||||
|
||||
GeneralMatrix& Q = dynareParams.getMatrixField(string("Sigma_e"));
|
||||
GeneralMatrix& H = dynareParams.getMatrixField(string("H"));
|
||||
#ifdef DEBUG
|
||||
Q.print();
|
||||
#endif
|
||||
GeneralMatrix& Hrtmp = dynareParams.getMatrixField(string("H"));
|
||||
GeneralMatrix * Hp;
|
||||
if (Hrtmp.numCols()==0 || Hrtmp.numRows()==0)
|
||||
{
|
||||
delete &Hrtmp;
|
||||
Hp = new GeneralMatrix(numVarobs,numVarobs);
|
||||
Hp->zeros();
|
||||
#ifdef DEBUG
|
||||
mexPrintf("finished local initialising of H \n");
|
||||
#endif
|
||||
}
|
||||
else
|
||||
Hp=&Hrtmp;
|
||||
|
||||
GeneralMatrix& H=*Hp;
|
||||
GeneralMatrix Y(data.numRows(),data.numCols());
|
||||
GeneralMatrix T(nr,nr);
|
||||
GeneralMatrix Z(numVarobs,nr);
|
||||
Z.zeros();
|
||||
for (int i = 0;i<numVarobs;++i)
|
||||
Z.get(i,mf[i]-1)=1;
|
||||
GeneralMatrix Pstar(nr,nr);
|
||||
GeneralMatrix R(nr,exo_nbr);
|
||||
GeneralMatrix ghx(endo_nbr,jcols-exo_nbr);
|
||||
GeneralMatrix ghx(endo_nbr,nsPred);
|
||||
GeneralMatrix ghu(endo_nbr,exo_nbr);
|
||||
|
||||
//Pinf=[]
|
||||
GeneralMatrix Pinf (nr,nr);
|
||||
Pinf.zeros();
|
||||
|
||||
double loglikelihood;
|
||||
try
|
||||
{
|
||||
|
||||
#ifdef DEBUG
|
||||
mexPrintf("Try construction of DsgeLikelihood\n");
|
||||
#endif
|
||||
DsgeLikelihood dl( a_init, Q, R,T, Z, Pstar, Pinf, H,data,Y,
|
||||
numPeriods, // const int INnumVarobs, // const int INnumTimeObs,
|
||||
order, endo_nbr, exo_nbr, nstatic, npred,nfwrd, num_of_observations,
|
||||
no_more_missing_observations, order_var, mfys, mf, xparam1,
|
||||
num_dp, deepParams, ub, lb, pshape, p6, p7, p3, p4, SteadyState, constant,
|
||||
dynareParams, dr, kstate, ghx, ghu, jcols, dfExt);
|
||||
dynareParams, dr, kstate, ghx, ghu, aux, iv, ic, jcols, dfExt);
|
||||
|
||||
double loglikelihood=dl.CalcLikelihood(xparam1);
|
||||
#ifdef DEBUG
|
||||
mexPrintf("Try CalcLikelihood\n");
|
||||
#endif
|
||||
loglikelihood=dl.CalcLikelihood(xparam1);
|
||||
|
||||
/*****************************************************************
|
||||
% OUTPUTS
|
||||
|
|
@ -161,6 +215,9 @@ mexFunction(int nlhs, mxArray *plhs[],
|
|||
% info : vector of informations about the penalty:
|
||||
% vll : vector of time-step log-likelihoods at xparam1.
|
||||
*****************************************************************/
|
||||
#ifdef DEBUG
|
||||
mexPrintf("Try Outputs with nper=%d, loglikelihood = %f\n",nper,loglikelihood);
|
||||
#endif
|
||||
if (nlhs >= 1)
|
||||
plhs[0] = mxCreateDoubleScalar(loglikelihood);
|
||||
if (nlhs >= 2)
|
||||
|
|
@ -169,23 +226,36 @@ mexFunction(int nlhs, mxArray *plhs[],
|
|||
{
|
||||
plhs[2] = mxCreateDoubleMatrix(endo_nbr, 1, mxREAL);
|
||||
Vector vss(mxGetPr(plhs[2]),endo_nbr);
|
||||
|
||||
#ifdef DEBUG
|
||||
mexPrintf("SteadyState size %d \n", dl.getSteadyState().length());
|
||||
dl.getSteadyState().print() ;
|
||||
mexPrintf("Try getSteadyState into vss size %d \n", vss.length());
|
||||
#endif
|
||||
vss= dl.getSteadyState();
|
||||
}
|
||||
/*********************
|
||||
if (nlhs >= 4)
|
||||
plhs[3] = mxCreateDoubleScalar((double)dl.getCostFlag());
|
||||
plhs[3] = mxCreateDoubleScalar((double)dl.getCostFlag());
|
||||
if (nlhs >= 5)
|
||||
plhs[4] = mxCreateDoubleMatrix(numVarobs,1, mxREAL);//dummy trend_coeff
|
||||
plhs[4] = mxCreateDoubleMatrix(numVarobs,1, mxREAL);//dummy trend_coeff
|
||||
if (nlhs >= 6)
|
||||
{
|
||||
// output full log-likelihood array
|
||||
/* Set the output pointer to the array of log likelihood. */
|
||||
// Set the output pointer to the array of log likelihood.
|
||||
std::vector<double>& vll=dl.getLikVector();
|
||||
plhs[5] = mxCreateDoubleMatrix(nper,1, mxREAL);
|
||||
double * mxll= mxGetPr(plhs[5]);
|
||||
// assign likelihood array
|
||||
for (int j=0;j<nper;++j)
|
||||
{
|
||||
mxll[j]=vll[j];
|
||||
#ifdef DEBUG
|
||||
mexPrintf("mxll[%d]=%f vll[%d]=%f\n",j, mxll[j], i, vll[j]);
|
||||
#endif
|
||||
}
|
||||
}
|
||||
*********************/
|
||||
}
|
||||
catch (const KordException &e)
|
||||
{
|
||||
|
|
|
|||
Loading…
Reference in New Issue