From d28442f2060ef9afe3d0bb237d57363981fec577 Mon Sep 17 00:00:00 2001 From: george Date: Fri, 20 Nov 2009 17:43:05 +0000 Subject: [PATCH] Update, still subject to further testing git-svn-id: https://www.dynare.org/svn/dynare/trunk@3160 ac1d8469-bf42-47a9-8791-bf33cf982152 --- mex/sources/estimation/DsgeLikelihood.cpp | 224 +++++++++++++----- mex/sources/estimation/DsgeLikelihood.h | 11 +- mex/sources/estimation/MexUtils.cpp | 24 +- mex/sources/estimation/MexUtils.h | 5 +- mex/sources/estimation/dr1.cpp | 3 + mex/sources/estimation/dynare_resolve.cpp | 130 ++++++---- .../estimation/tests/DsgeLikelihood_mexf.cpp | 112 +++++++-- 7 files changed, 378 insertions(+), 131 deletions(-) diff --git a/mex/sources/estimation/DsgeLikelihood.cpp b/mex/sources/estimation/DsgeLikelihood.cpp index a5ab94a4e..3b44265c2 100644 --- a/mex/sources/estimation/DsgeLikelihood.cpp +++ b/mex/sources/estimation/DsgeLikelihood.cpp @@ -36,7 +36,8 @@ DsgeLikelihood::DsgeLikelihood( Vector& inA_init, GeneralMatrix& inQ, GeneralMa const Vector& INub, const Vector& INlb, const vector&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&INiv, vector&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* vll=new std::vector (nper); - vll=new std::vector (numTimeObs);// vector of likelihoods +// vll=new std::vector (numTimeObs);// vector of likelihoods + vll=new vector(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 *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&iv= dynareParams.getIntVectorField(string("restrict_var_list")); - vector&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&lgyidx2varobs= dynareParams.getIntVectorField(string("lgyidx2varobs")); for (i=0;i&iv; //= dynareParams.getIntVectorField(string("restrict_var_list")); + vector⁣ //= 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&iv,vector&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&iv, vector&ic ,const int jcols, const char *dfExt); //, KordpDynare& inModel, Approximation& INapprox ); DsgeLikelihood( const Vector¶ms,const GeneralMatrix&data, const vector& data_index, const int gend, const int number_of_observations, const bool no_more_missing_observations);//, KordpDynare& model ); // constructor, and diff --git a/mex/sources/estimation/MexUtils.cpp b/mex/sources/estimation/MexUtils.cpp index 2d4318c45..aaaa12ec3 100644 --- a/mex/sources/estimation/MexUtils.cpp +++ b/mex/sources/estimation/MexUtils.cpp @@ -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::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 } diff --git a/mex/sources/estimation/MexUtils.h b/mex/sources/estimation/MexUtils.h index 51de463d9..7a1ccdf37 100644 --- a/mex/sources/estimation/MexUtils.h +++ b/mex/sources/estimation/MexUtils.h @@ -73,7 +73,10 @@ class MexStruct :public virtual GeneralParams { map::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() ); } diff --git a/mex/sources/estimation/dr1.cpp b/mex/sources/estimation/dr1.cpp index 7589486d1..dda768cad 100644 --- a/mex/sources/estimation/dr1.cpp +++ b/mex/sources/estimation/dr1.cpp @@ -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 } } diff --git a/mex/sources/estimation/dynare_resolve.cpp b/mex/sources/estimation/dynare_resolve.cpp index 1a21ecc89..c9ce2e951 100644 --- a/mex/sources/estimation/dynare_resolve.cpp +++ b/mex/sources/estimation/dynare_resolve.cpp @@ -59,24 +59,38 @@ DsgeLikelihood::dynareResolveDR(vector&iv,vector&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&iv,vector&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&iv,vector&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(); - - vectori_n_iv(n_iv); - for (int i=0;ii_n_iv(n_iv); + for (int i=0;i 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;ispan2end(A.numRows()-n_iv); - for (int i=0;ispan2end(A.numRows()-n_iv); + for (int i=0;i&mfys=dynareParams.getIntVectorField("mfys"); vector&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&order_var = dr.getIntVectorField(string("order_var")); +#ifdef DEBUG + for ( gg=0;gg&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&iv= dynareParams.getIntVectorField(string("restrict_var_list")); + vector&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= 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& vll=dl.getLikVector(); plhs[5] = mxCreateDoubleMatrix(nper,1, mxREAL); double * mxll= mxGetPr(plhs[5]); // assign likelihood array for (int j=0;j