Update: jacobian reordering methods and update of sparse tensor population. Also, test .cpp using new 1 lag 1 lead model derived from fs2000a.
git-svn-id: https://www.dynare.org/svn/dynare/trunk@2356 ac1d8469-bf42-47a9-8791-bf33cf982152time-shift
george
parent
a750cba621
commit
a39c2eb17f
|
|
@ -13,6 +13,7 @@
|
|||
*/
|
||||
// GP, based on work by O.Kamenik
|
||||
|
||||
#include <vector>
|
||||
#include "first_order.h"
|
||||
#include "k_ord_dynare.h"
|
||||
|
||||
|
|
@ -79,11 +80,12 @@ 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)
|
||||
Journal& jr, DynamicModelDLL& dynamicDLL, double sstol, const int* var_order,
|
||||
const TwoDMatrix * llincidence )
|
||||
: 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)
|
||||
md(1), dnl(NULL), denl(NULL), dsnl(NULL), ss_tol(sstol), varOrder( var_order), ll_Incidence(llincidence)
|
||||
{
|
||||
#ifdef DEBUG
|
||||
mexPrintf("k_ord_dynare Dynare constructor: ny=%d, order=%d, nPar=%d .\n", nY,nOrder,nPar);
|
||||
|
|
@ -109,6 +111,8 @@ KordpDynare::KordpDynare(const char** endo, int num_endo,
|
|||
vCov = new TwoDMatrix(*(vCov));
|
||||
******/
|
||||
// throw DynareException(__FILE__, __LINE__, string("Could not open model file ")+modName);
|
||||
JacobianIndices= ReorderDynareJacobianIndices( varOrder);
|
||||
|
||||
}
|
||||
catch (...){
|
||||
mexPrintf("k_ord_dynare: dynare constructor, error in StateNamelist construction.\n");
|
||||
|
|
@ -126,7 +130,9 @@ KordpDynare::KordpDynare(const KordpDynare& dynare)
|
|||
nSteps(dynare.nSteps), nOrder(dynare.nOrder), journal(dynare.journal),
|
||||
dynamicDLL(dynare.dynamicDLL), //modName(dynare.modName),
|
||||
ySteady(NULL), params(NULL), vCov(NULL), md(dynare.md),
|
||||
dnl(NULL), denl(NULL), dsnl(NULL), ss_tol(dynare.ss_tol)
|
||||
dnl(NULL), denl(NULL), dsnl(NULL), ss_tol(dynare.ss_tol),
|
||||
varOrder(dynare.varOrder), ll_Incidence(dynare.ll_Incidence),
|
||||
JacobianIndices(dynare.JacobianIndices)
|
||||
{
|
||||
/// model = dynare.model->clone();
|
||||
ySteady = new Vector(*(dynare.ySteady));
|
||||
|
|
@ -236,12 +242,12 @@ void KordpDynare::calcDerivatives(const Vector& yy, const Vector& xx)
|
|||
mexPrintf("k_ord_dynare.cpp: Error in dynamicDLL.eval in calcDerivatives");
|
||||
return;
|
||||
}
|
||||
|
||||
if ((nJcols!=g1->ncols()) && ( nY != g1->nrows())) {
|
||||
mexPrintf("k_ord_dynare.cpp: Error in calcDerivatives: dynamicDLL.eval returned wrong jacobian");
|
||||
return;
|
||||
}
|
||||
|
||||
ReorderCols(g1, JacobianIndices);
|
||||
// ReorderCols(out, varOrder);
|
||||
#ifdef DEBUG
|
||||
mexPrintf("k_ord_dynare.cpp: populate FSSparseTensor in calcDerivatives: cols=%d , rows=%d\n"
|
||||
, g1->ncols(),g1->nrows());
|
||||
|
|
@ -249,17 +255,21 @@ void KordpDynare::calcDerivatives(const Vector& yy, const Vector& xx)
|
|||
|
||||
// model derivatives FSSparseTensor instance for single order only
|
||||
//(higher orders requires Symetry to insert in particular position.)
|
||||
FSSparseTensor *mdTi=(new FSSparseTensor (1, g1->ncols(),g1->nrows()));
|
||||
FSSparseTensor *mdTi=(new FSSparseTensor (1, g1->ncols(),g1->nrows()));
|
||||
|
||||
for (int i = 0; i<g1->ncols(); i++){
|
||||
for (int j = 0; j<g1->nrows(); j++){
|
||||
if (g1->get(j,i)!=0.0) // populate sparse if not zero
|
||||
mdTi->insert(i, j,g1->get(j,i));
|
||||
if (g1->get(j,i)!=0.0){ // populate sparse if not zero
|
||||
IntSequence s(1,0);
|
||||
s[0]=i;
|
||||
mdTi->insert(s, j,g1->get(j,i));
|
||||
}
|
||||
}
|
||||
}
|
||||
// md container
|
||||
// md=*(new TensorContainer<FSSparseTensor>(1));
|
||||
// FSSparseTensor mdSTi(mdTi);
|
||||
md.clear();
|
||||
md.clear();// this is to be used only for 1st order!!
|
||||
md.insert(mdTi);//(&mdTi);
|
||||
}
|
||||
|
||||
|
|
@ -350,6 +360,108 @@ void KordpDynare::writeModelInfo(Journal& jr) const
|
|||
}
|
||||
|
||||
|
||||
/************************************
|
||||
* Reorder DynareJacobianIndices of variables in a vector according to
|
||||
* given int * varOrder together with lead & lag incidence matrix and
|
||||
* any the extra columns for exogenous vars, and then,
|
||||
* reorders its blocks given by the varOrder and the Dynare++ expectations:
|
||||
|
||||
* extra nboth+ npred (t-1) lags
|
||||
* varOrder
|
||||
static:
|
||||
pred
|
||||
both
|
||||
forward
|
||||
* extra both + nforw (t+1) leads, and
|
||||
* extra exogen
|
||||
|
||||
* so to match the jacobian organisation expected by the Appoximation class
|
||||
both + nforw (t+1) leads
|
||||
static
|
||||
pred
|
||||
both
|
||||
forward
|
||||
nboth+ npred (t-1) lags
|
||||
exogen
|
||||
************************************/
|
||||
|
||||
int * KordpDynare::ReorderDynareJacobianIndices( const int * 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));
|
||||
vector <int> tmp(nY);
|
||||
int i,j, rjoff=nJcols-nExog-1; //, ll_off, j;
|
||||
try{
|
||||
for (int ll_row=0; ll_row< ll_Incidence->nrows(); ll_row++)
|
||||
{
|
||||
// reorder in orde-var order & populate temporary nEndo (sparse) vector with
|
||||
// the lag, current and lead blocks of the jacobian respectively
|
||||
for (i=0;i<nY;i++){
|
||||
// tmp[varOrder[j]]=(int)ll_Incidence->get(ll_row,j);
|
||||
tmp[i]=((int)ll_Incidence->get(ll_row,varOrder[i]-1));
|
||||
}
|
||||
// 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;
|
||||
rjoff--;
|
||||
if (rjoff<0){
|
||||
// mexPrintf(" Error in ReorderIndices - negative rjoff index?");
|
||||
break;
|
||||
// return NULL;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
} catch (const TLException& e) {
|
||||
mexPrintf("Caugth TL exception in ReorderIndices: ");
|
||||
e.print();
|
||||
return NULL;// 255;
|
||||
}catch (...){
|
||||
mexPrintf(" Error in ReorderIndices - wrong index?");
|
||||
}
|
||||
//add the indices for the nExog exogenous jacobians
|
||||
for (j=nJcols-nExog;j<nJcols;j++){
|
||||
JacobianIndices[j]=j;
|
||||
}
|
||||
#ifdef DEBUG
|
||||
for (j=0;j<nJcols;j++)
|
||||
mexPrintf("ReorderDynareJacobianIndice: [%d] =%d .\n", j, JacobianIndices[j]);
|
||||
#endif
|
||||
return JacobianIndices;
|
||||
}
|
||||
|
||||
|
||||
/************************************
|
||||
* Reorder first set of columns of variables in a (jacobian) matrix
|
||||
* according to order given in varsOrder together with the extras
|
||||
* assuming tdx ncols() - nExog is eaqual or less than length of varOrder and
|
||||
* of any of its elements too.
|
||||
************************************/
|
||||
|
||||
void KordpDynare::ReorderCols(TwoDMatrix * tdx, const int * varsOrder){
|
||||
|
||||
TwoDMatrix tmp(*tdx); // temporary 2D matrix
|
||||
TwoDMatrix &tmpR=tmp;
|
||||
tdx->zeros();// empty original matrix
|
||||
// reorder the columns
|
||||
try{
|
||||
for (int i =0; i<tdx->ncols() ; i++)
|
||||
tdx->copyColumn(tmpR,varsOrder[i],i);
|
||||
} catch (const TLException& e) {
|
||||
printf("Caugth TL exception in ReorderColumns: ");
|
||||
e.print();
|
||||
return;// 255;
|
||||
}catch (...){
|
||||
mexPrintf(" Error in ReorderColumns - wrong index?");
|
||||
}
|
||||
}
|
||||
|
||||
/****************************
|
||||
* K-Order Perturbation instance of Jacobian:
|
||||
************************************/
|
||||
|
|
|
|||
|
|
@ -135,12 +135,16 @@ friend class DynamicModelDLL;
|
|||
DynareExogNameList* denl;
|
||||
DynareStateNameList* dsnl;
|
||||
const double ss_tol;
|
||||
const int* varOrder;
|
||||
const TwoDMatrix * ll_Incidence;
|
||||
int * 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);
|
||||
Journal& jr, DynamicModelDLL& dynamicDLL, double sstol,
|
||||
const int* varOrder, const TwoDMatrix * ll_Incidence );
|
||||
|
||||
/** Makes a deep copy of the object. */
|
||||
KordpDynare(const KordpDynare& dyn);
|
||||
|
|
@ -201,8 +205,11 @@ public:
|
|||
/// void writeDump(const std::string& basename) const;
|
||||
DynamicModel* clone() const
|
||||
{return new KordpDynare(*this);}
|
||||
void ReorderCols(TwoDMatrix * tdx, const int * varOrder);
|
||||
private:
|
||||
void writeModelInfo(Journal& jr) const;
|
||||
int * ReorderDynareJacobianIndices( const int * varOrder);
|
||||
void ReorderBlocks(TwoDMatrix * tdx, const int * varOrder);
|
||||
};
|
||||
|
||||
/****************************
|
||||
|
|
|
|||
|
|
@ -144,7 +144,7 @@ extern "C" {
|
|||
mxFldp = mxGetField(dr, 0,"nstatic" );
|
||||
const int nStat = (int)mxGetScalar(mxFldp);
|
||||
mxFldp = mxGetField(dr, 0,"npred" );
|
||||
const int nPred = (int)mxGetScalar(mxFldp);
|
||||
int nPred = (int)mxGetScalar(mxFldp);
|
||||
mxFldp = mxGetField(dr, 0,"nspred" );
|
||||
const int nsPred = (int)mxGetScalar(mxFldp);
|
||||
mxFldp = mxGetField(dr, 0,"nboth" );
|
||||
|
|
@ -164,6 +164,22 @@ extern "C" {
|
|||
mxFldp = mxGetField(M_, 0,"maximum_lag" );
|
||||
const int nMax_lag = (int)mxGetScalar(mxFldp);
|
||||
|
||||
nPred -= nBoth; // correct nPred for nBoth.
|
||||
|
||||
mxFldp = mxGetField(dr, 0,"order_var" );
|
||||
int * var_order = (int *) 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);
|
||||
|
||||
mxFldp = mxGetField(M_, 0,"lead_lag_incidence" );
|
||||
dparams = (double *) mxGetData(mxFldp);
|
||||
npar = (int)mxGetN(mxFldp);
|
||||
TwoDMatrix * llincidence = new TwoDMatrix(npar, nEndo, dparams);
|
||||
|
||||
const int jcols = nExog+nEndo+nsPred+nsForw; // Num of Jacobian columns
|
||||
mexPrintf("k_order_perturbation: jcols= %d .\n", jcols);
|
||||
|
||||
|
|
@ -236,7 +252,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);
|
||||
jcols, nSteps, kOrder, journal, dynamicDLL, sstol, var_order,
|
||||
llincidence );
|
||||
/************
|
||||
// make list of shocks for which we will do IRFs
|
||||
vector<int> irf_list_ind;
|
||||
|
|
|
|||
|
|
@ -76,7 +76,8 @@ int main(int argc, char* argv[])
|
|||
// };
|
||||
const int nSteady = 31;//29, 16 (int)mxGetM(mxFldp);
|
||||
Vector * ySteady = new Vector(dYSparams, nSteady);
|
||||
|
||||
|
||||
|
||||
//mxFldp = mxGetField(dr, 0,"nstatic" );
|
||||
const int nStat = 7;//(int)mxGetScalar(mxFldp);
|
||||
// mxFldp = mxGetField(dr, 0,"npred" );
|
||||
|
|
@ -99,15 +100,52 @@ int main(int argc, char* argv[])
|
|||
// it_ should be set to M_.maximum_lag
|
||||
//mxFldp = mxGetField(M_, 0,"maximum_lag" );
|
||||
const int nMax_lag = 1;//(int)mxGetScalar(mxFldp);
|
||||
|
||||
const int jcols = nExog+nEndo+nsPred+nsForw; // Num of Jacobian columns
|
||||
|
||||
int var_order[]//[18]
|
||||
= {
|
||||
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);
|
||||
|
||||
const double ll_incidence []//[3][18]
|
||||
= {
|
||||
/* 1, 2, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 4, 5
|
||||
, 6, 0, 0, 0
|
||||
, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20
|
||||
, 21, 22, 23, 24
|
||||
, 25, 26, 27, 28, 0, 0, 0, 0, 29, 0, 0, 0, 0, 0
|
||||
, 0, 0, 30, 31
|
||||
*/
|
||||
1, 7, 25
|
||||
, 2, 8, 26
|
||||
, 0, 9, 27
|
||||
, 0, 10, 28
|
||||
, 0, 11, 0
|
||||
, 0, 12, 0
|
||||
, 3, 13, 0
|
||||
, 0, 14, 0
|
||||
, 0, 15, 29
|
||||
, 0, 16, 0
|
||||
, 0, 17, 0
|
||||
, 0, 18, 0
|
||||
, 4, 19, 0
|
||||
, 5, 20, 0
|
||||
, 6, 21, 0
|
||||
, 0, 22, 0
|
||||
, 0, 23, 30
|
||||
, 0, 24, 31
|
||||
};
|
||||
TwoDMatrix * llincidence = new TwoDMatrix ( 3, nEndo, ll_incidence );
|
||||
|
||||
const int jcols = nExog+nEndo+nsPred+nsForw; // Num of Jacobian columns
|
||||
#ifdef DEBUG
|
||||
mexPrintf("k_order_perturbation: jcols= %d .\n", jcols);
|
||||
|
||||
#endif
|
||||
//mxFldp= mxGetField(M_, 0,"endo_names" );
|
||||
const int nendo = 18;//16(int)mxGetM(mxFldp);
|
||||
const int widthEndo = 6;// (int)mxGetN(mxFldp);
|
||||
const char * cNamesCharStr= "mPceWRkdnlggYPydPc yp__ A22 __oo oobb bbss ss ";
|
||||
// const char** endoNamesMX= DynareMxArrayToString( mxFldp,nendo,widthEndo);
|
||||
// const char** endoNamesMX= DynareMxArrayToString( mxFldp,nendo,widthEndo);
|
||||
const char** endoNamesMX= DynareMxArrayToString( cNamesCharStr, nendo, widthEndo);
|
||||
#ifdef DEBUG
|
||||
for (int i = 0; i < nEndo; i++) {
|
||||
|
|
@ -120,7 +158,6 @@ int main(int argc, char* argv[])
|
|||
// const char** exoNamesMX= DynareMxArrayToString( mxFldp,nexo,widthExog);
|
||||
const char * cExoNamesCharStr= "ee__am";
|
||||
const char** exoNamesMX= DynareMxArrayToString( cExoNamesCharStr,nexo,widthExog);
|
||||
|
||||
#ifdef DEBUG
|
||||
for (int i = 0; i < nexo; i++) {
|
||||
mexPrintf("k_ord_perturbation: ExoNameList[%d][0]= %s.\n", i, exoNamesMX[i] );
|
||||
|
|
@ -135,7 +172,7 @@ int main(int argc, char* argv[])
|
|||
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] ); }
|
||||
//, , 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 < nSteady; i++) {
|
||||
|
|
@ -149,7 +186,7 @@ int main(int argc, char* argv[])
|
|||
const int nSteps =0; // Dynare++ solving steps, for time being default to 0 = deterministic steady state
|
||||
const double sstol=1.e-13; //NL solver tolerance from
|
||||
|
||||
THREAD_GROUP::max_parallel_threads = 2;//params.num_threads;
|
||||
THREAD_GROUP::max_parallel_threads = 1;//2 params.num_threads;
|
||||
|
||||
try {
|
||||
// make journal name and journal
|
||||
|
|
@ -165,9 +202,10 @@ int main(int argc, char* argv[])
|
|||
mexPrintf("k_order_perturbation: Calling dynare constructor.\n");
|
||||
#endif
|
||||
// make KordpDynare object
|
||||
KordpDynare dynare(endoNamesMX, nEndo, exoNamesMX, nExog, nPar, // paramNames,
|
||||
ySteady, vCov, modParams, nStat, nPred, nForw, nBoth,
|
||||
jcols, nSteps, kOrder, journal, dynamicDLL, sstol);
|
||||
KordpDynare dynare(endoNamesMX, nEndo, exoNamesMX, nExog, nPar // paramNames,
|
||||
, ySteady, vCov, modParams, nStat, nPred, nForw, nBoth
|
||||
, jcols, nSteps, kOrder, journal, dynamicDLL, sstol, var_order
|
||||
, llincidence );
|
||||
try {
|
||||
// intiate tensor library
|
||||
#ifdef DEBUG
|
||||
|
|
@ -176,7 +214,6 @@ int main(int argc, char* argv[])
|
|||
tls.init(dynare.order(),
|
||||
dynare.nstat()+2*dynare.npred()+3*dynare.nboth()+
|
||||
2*dynare.nforw()+dynare.nexog());
|
||||
|
||||
// construct main K-order approximation class
|
||||
// FistOrderApproximation app(dynare, journal, nSteps);
|
||||
#ifdef DEBUG
|
||||
|
|
@ -202,11 +239,9 @@ int main(int argc, char* argv[])
|
|||
mexErrMsgTxt("Caught Sylv exception.");
|
||||
}
|
||||
|
||||
|
||||
// get latest ysteady
|
||||
double * dYsteady = (dynare.getSteady().base());
|
||||
ySteady = (Vector*)(&dynare.getSteady());
|
||||
|
||||
} catch (const KordException& e) {
|
||||
printf("Caugth Kord exception: ");
|
||||
e.print();
|
||||
|
|
|
|||
Loading…
Reference in New Issue