diff --git a/mex/sources/korderpert/src/k_ord_dynare.cpp b/mex/sources/korderpert/src/k_ord_dynare.cpp index 0fd080364..9672888c8 100644 --- a/mex/sources/korderpert/src/k_ord_dynare.cpp +++ b/mex/sources/korderpert/src/k_ord_dynare.cpp @@ -86,12 +86,13 @@ KordpDynare::KordpDynare(const char** endo, int num_endo, const char** exo, int nexog, int npar, //const char** par, Vector* ysteady, TwoDMatrix* vcov, Vector* inParams, int nstat, int npred, int nforw, int nboth, const int jcols, const int nsteps, int norder, //const char* modName, - Journal& jr, DynamicModelDLL& dynamicDLL, double sstol, const int* var_order, - const TwoDMatrix * llincidence ) + Journal& jr, DynamicModelDLL& dynamicDLL, double sstol, const vector* var_order, + const TwoDMatrix * llincidence, double criterium ) : nStat(nstat), nBoth(nboth), nPred(npred), nForw(nforw), nExog(nexog), nPar(npar), nYs(npred + nboth), nYss(nboth + nforw),nY(num_endo), nJcols(jcols), nSteps(nsteps), nOrder(norder), journal(jr), dynamicDLL(dynamicDLL), ySteady(ysteady), vCov(vcov), params (inParams), - md(1), dnl(NULL), denl(NULL), dsnl(NULL), ss_tol(sstol), varOrder( var_order), ll_Incidence(llincidence) + md(1), dnl(NULL), denl(NULL), dsnl(NULL), ss_tol(sstol), varOrder( var_order), + ll_Incidence(llincidence), qz_criterium(criterium) { #ifdef DEBUG mexPrintf("k_ord_dynare Dynare constructor: ny=%d, order=%d, nPar=%d .\n", nY,nOrder,nPar); @@ -138,7 +139,7 @@ KordpDynare::KordpDynare(const KordpDynare& dynare) ySteady(NULL), params(NULL), vCov(NULL), md(dynare.md), dnl(NULL), denl(NULL), dsnl(NULL), ss_tol(dynare.ss_tol), varOrder(dynare.varOrder), ll_Incidence(dynare.ll_Incidence), - JacobianIndices(dynare.JacobianIndices) + JacobianIndices(dynare.JacobianIndices), qz_criterium(dynare.qz_criterium) { /// model = dynare.model->clone(); ySteady = new Vector(*(dynare.ySteady)); @@ -391,7 +392,7 @@ Vector * KordpDynare::LLxSteady( const Vector& yS){ // populate (non-sparse) vector with ysteady values for (int i=0;iget(ll_row,i)) - (*llxSteady)[ll_Incidence->get(ll_row,i)-1] = yS[i]; + (*llxSteady)[((int)ll_Incidence->get(ll_row,i))-1] = yS[i]; } } } catch (const TLException& e) { @@ -435,16 +436,20 @@ Vector * KordpDynare::LLxSteady( const Vector& yS){ exogen ************************************/ -int * KordpDynare::ReorderDynareJacobianIndices( const int * varOrder){ +vector * KordpDynare::ReorderDynareJacobianIndices( const vector * varOrder){ // if ((nJcols != tdx->ncols()) && ( nY != tdx->nrows())) { // mexPrintf("k_ord_dynare.cpp: Error in ReorderBlocks: wrong size of jacobian"); // return; // } // create temporary square 2D matrix size nEndo x nEndo (sparse) // for the lag, current and lead blocks of the jacobian - int * JacobianIndices = (int*) calloc(nJcols+1, sizeof(int)); +// int * JacobianIndices = (int*) calloc(nJcols+1, sizeof(int)); + vector * JacobianIndices = new vector(nJcols); vector tmp(nY); int i,j, rjoff=nJcols-nExog-1; //, ll_off, j; +#ifdef DEBUG + mexPrintf("ReorderDynareJacobianIndice:ll_Incidence->nrows() =%d .\n", ll_Incidence->nrows()); +#endif try{ for (int ll_row=0; ll_row< ll_Incidence->nrows(); ll_row++) { @@ -452,13 +457,17 @@ int * KordpDynare::ReorderDynareJacobianIndices( const int * varOrder){ // the lag, current and lead blocks of the jacobian respectively for (i=0;iget(ll_row,j); - tmp[i]=((int)ll_Incidence->get(ll_row,varOrder[i]-1)); + tmp[i]=((int)ll_Incidence->get(ll_row,(*varOrder)[i]-1)); +#ifdef DEBUG + mexPrintf("get(ll_row,(*varOrder)[%d]-1)) = tmp[%d]=%d .\n", + i, i, (int)ll_Incidence->get(ll_row,(*varOrder)[i]-1)); +#endif } // write the reordered blocks back to the jacobian // in reverse order for (j=nY-1;j>=0;j--){ if (tmp[j]){ - JacobianIndices[rjoff]=tmp[j] -1; + (*JacobianIndices)[rjoff]=tmp[j] -1; rjoff--; if (rjoff<0){ // mexPrintf(" Error in ReorderIndices - negative rjoff index?"); @@ -477,11 +486,11 @@ int * KordpDynare::ReorderDynareJacobianIndices( const int * varOrder){ } //add the indices for the nExog exogenous jacobians for (j=nJcols-nExog;j * vOrder){ + + if (tdx->ncols() > vOrder->size()){ + mexPrintf(" Error in ReorderColumns - size of order var is too small"); + return; + } + TwoDMatrix tmp(*tdx); // temporary 2D matrix + TwoDMatrix &tmpR=tmp; + tdx->zeros();// empty original matrix + // reorder the columns + try{ + for (int i =0; incols() ; i++) + tdx->copyColumn(tmpR,(*vOrder)[i],i); + } catch (const TLException& e) { + printf("Caugth TL exception in ReorderColumns: "); + e.print(); + return;// 255; + }catch (...){ + mexPrintf(" Error in ReorderColumns - wrong index?"); + } +} +void KordpDynare::ReorderCols(TwoDMatrix * tdx, const int * vOrder){ TwoDMatrix tmp(*tdx); // temporary 2D matrix TwoDMatrix &tmpR=tmp; @@ -502,7 +532,7 @@ void KordpDynare::ReorderCols(TwoDMatrix * tdx, const int * varsOrder){ // reorder the columns try{ for (int i =0; incols() ; i++) - tdx->copyColumn(tmpR,varsOrder[i],i); + tdx->copyColumn(tmpR,vOrder[i],i); } catch (const TLException& e) { printf("Caugth TL exception in ReorderColumns: "); e.print(); diff --git a/mex/sources/korderpert/src/k_ord_dynare.h b/mex/sources/korderpert/src/k_ord_dynare.h index 6c32cb23e..fc962ffdb 100644 --- a/mex/sources/korderpert/src/k_ord_dynare.h +++ b/mex/sources/korderpert/src/k_ord_dynare.h @@ -141,16 +141,18 @@ friend class DynamicModelDLL; DynareExogNameList* denl; DynareStateNameList* dsnl; const double ss_tol; - const int* varOrder; + const vector* varOrder; const TwoDMatrix * ll_Incidence; - int * JacobianIndices; + double qz_criterium; + vector * JacobianIndices; public: KordpDynare(const char** endo, int num_endo, const char** exo, int num_exo, int num_par, //const char** par, Vector* ySteady, TwoDMatrix* vCov, Vector* params, int nstat,int nPred, int nforw, int nboth, const int nJcols, const int nSteps, const int ord, //const char* modName, Journal& jr, DynamicModelDLL& dynamicDLL, double sstol, - const int* varOrder, const TwoDMatrix * ll_Incidence ); + const vector* varOrder, const TwoDMatrix * ll_Incidence, + double qz_criterium ); /** Makes a deep copy of the object. */ KordpDynare(const KordpDynare& dyn); @@ -212,11 +214,13 @@ public: DynamicModel* clone() const {return new KordpDynare(*this);} void ReorderCols(TwoDMatrix * tdx, const int * varOrder); + void ReorderCols(TwoDMatrix * tdx, const vector * varOrder); Vector * KordpDynare::LLxSteady( const Vector& yS); // returns ySteady extended with leads and lags private: void writeModelInfo(Journal& jr) const; int * ReorderDynareJacobianIndices( const int * varOrder); + vector * ReorderDynareJacobianIndices( const vector * varOrder); void ReorderBlocks(TwoDMatrix * tdx, const int * varOrder); }; diff --git a/mex/sources/korderpert/src/k_order_perturbation.cpp b/mex/sources/korderpert/src/k_order_perturbation.cpp index f68dfc080..2f05c7373 100644 --- a/mex/sources/korderpert/src/k_order_perturbation.cpp +++ b/mex/sources/korderpert/src/k_order_perturbation.cpp @@ -110,6 +110,11 @@ extern "C" { kOrder = (int)mxGetScalar(mxFldp); else kOrder = 1; + + double qz_criterium = 1+1e-6; + mxFldp = mxGetField(options_, 0,"qz_criterium" ); + if (mxIsNumeric(mxFldp)) + qz_criterium = (int)mxGetScalar(mxFldp); mxFldp = mxGetField(M_, 0,"params" ); double * dparams = (double *) mxGetData(mxFldp); @@ -117,30 +122,25 @@ extern "C" { Vector * modParams = new Vector(dparams, npar); #ifdef DEBUG mexPrintf("k_ord_perturbation: nParams=%d .\n",npar); - for (int i = 0; i < npar; i++) { - mexPrintf("k_ord_perturbation: dParams[%d]= %g.\n", i, dparams+i*(sizeof(double)) ); } for (int i = 0; i < npar; i++) { mexPrintf("k_ord_perturbation: Params[%d]= %g.\n", i, (*modParams)[i]); } - #endif - const mxArray* const mxParFldp = mxGetField(M_, 0,"params" ); - Vector params_vec((const double*)mxGetPr(mxParFldp), npar); + //const mxArray* const mxParFldp = mxGetField(M_, 0,"params" ); + //Vector params_vec((const double*)mxGetPr(mxParFldp), npar); #ifdef DEBUG - for (int i = 0; i < npar; i++) { - mexPrintf("k_ord_perturbation: params_vec[%d]= %g.\n", i, params_vec[i] ); } - +// for (int i = 0; i < npar; i++) { +// mexPrintf("k_ord_perturbation: params_vec[%d]= %g.\n", i, params_vec[i] ); } #endif - mxFldp = mxGetField(M_, 0,"Sigma_e" ); dparams = (double *) mxGetData(mxFldp); npar = (int)mxGetN(mxFldp); TwoDMatrix * vCov = new TwoDMatrix(npar, npar, dparams); -// mxFldp = mxGetField(oo_, 0,"steady_state" ); // use in order of declaration + mxFldp = mxGetField(oo_, 0,"steady_state" ); // use in order of declaration // mxFldp = mxGetField(dr, 0,"ys" ); // and not in order of dr.order_var - mxFldp = mxGetField(oo_, 0,"dyn_ys" ); // extended ys +// mxFldp = mxGetField(oo_, 0,"dyn_ys" ); // and NOT extended ys dparams = (double *) mxGetData(mxFldp); const int nSteady = (int)mxGetM(mxFldp); Vector * ySteady = new Vector(dparams, nSteady); @@ -172,18 +172,39 @@ extern "C" { nPred -= nBoth; // correct nPred for nBoth. mxFldp = mxGetField(dr, 0,"order_var" ); - int * var_order = (int *) mxGetData(mxFldp); + dparams = (double *) mxGetData(mxFldp); npar = (int)mxGetM(mxFldp); if (npar != nEndo) { //(nPar != npar) mexErrMsgTxt("Incorrect number of input var_order vars.\n"); //return; } -// Vector * varOrder = new Vector(var_order, nEndo); + vector * var_order_vp = (new vector(nEndo));//nEndo)); + for (int v=0;vget(%d,%d) =%d .\n", + j, i, (int)llincidence->get(j,i));}} +#endif const int jcols = nExog+nEndo+nsPred+nsForw; // Num of Jacobian columns mexPrintf("k_order_perturbation: jcols= %d .\n", jcols); @@ -218,16 +239,14 @@ extern "C" { ************/ if ((nEndo != nendo)||(nExog != nexo)) { //(nPar != npar) mexErrMsgTxt("Incorrect number of input parameters.\n"); - //return; + return; } #ifdef DEBUG for (int i = 0; i < nEndo; i++) { mexPrintf("k_ord_perturbation: EndoNameList[%d]= %s.\n", i, endoNamesMX[i] ); } - for (int i = 0; i < nPar; i++) { - mexPrintf("k_ord_perturbation: params_vec[%d]= %g.\n", i, params_vec[i] ); } - for (int i = 0; i < nPar; i++) { - mexPrintf("k_ord_perturbation: Params[%d]= %g.\n", i, (*modParams)[i]); } +// for (int i = 0; i < nPar; i++) { +// mexPrintf("k_ord_perturbation: Params[%d]= %g.\n", i, (*modParams)[i]); } for (int i = 0; i < nSteady; i++) { mexPrintf("k_ord_perturbation: ysteady[%d]= %g.\n", i, (*ySteady)[i]); } @@ -257,8 +276,8 @@ extern "C" { // make KordpDynare object KordpDynare dynare(endoNamesMX, nEndo, exoNamesMX, nExog, nPar, // paramNames, ySteady, vCov, modParams, nStat, nPred, nForw, nBoth, - jcols, nSteps, kOrder, journal, dynamicDLL, sstol, var_order, - llincidence ); + jcols, nSteps, kOrder, journal, dynamicDLL, sstol, var_order_vp, + llincidence,qz_criterium ); /************ // make list of shocks for which we will do IRFs vector irf_list_ind; @@ -650,8 +669,8 @@ void DynamicModelDLL::eval(const Vector&y, const TwoDMatrix&x, const Vector* mo , g1->ncols(),g1->nrows()); for (int i = 0; i < modParams->length(); i++) { mexPrintf("k_ord_perturbation: Params[%d]= %g.\n", i, (*modParams)[i]); } - for (int i = 0; i < length; i++) { - mexPrintf("k_ord_perturbation: Params[%d]= %g.\n", i, y[i]);} + for (int i = 0; i < jcols-nExog; i++) { + mexPrintf("k_ord_perturbation: Ys[%d]= %g.\n", i, y[i]);} mexPrintf("k_order_perturbation: call Dynamic dParams= %d , , dy = %d dx = %d .\n" ,dbParams[0],dy[0],dx[0]); diff --git a/mex/sources/korderpert/src/k_order_test_main.cpp b/mex/sources/korderpert/src/k_order_test_main.cpp index f8c9487cf..4eec39308 100644 --- a/mex/sources/korderpert/src/k_order_test_main.cpp +++ b/mex/sources/korderpert/src/k_order_test_main.cpp @@ -31,7 +31,7 @@ int main(int argc, char* argv[]) { - + double qz_criterium = 1+1e-6; const int check_flag = 0; const char* fName = "fs2000k";//mxArrayToString(mFname); @@ -48,7 +48,6 @@ int main(int argc, char* argv[]) 0.7870, 0.0200 }; - Vector * modParams = new Vector(dparams, npar); #ifdef DEBUG @@ -106,7 +105,10 @@ int main(int argc, char* argv[]) // 5, 6, 8, 10, 11, 12, 16, 7, 13, 14, 15, 1, 2, 3, 4, 9, 17, 18 }; //Vector * varOrder = new Vector(var_order, nEndo); - + vector * var_order_vp = new vector(nEndo);//nEndo)); + for (int v=0;v