diff --git a/mex/sources/korderpert/src/fs2000k.mod b/mex/sources/korderpert/src/fs2000k.mod index 072bb231f..9ba2b6061 100644 --- a/mex/sources/korderpert/src/fs2000k.mod +++ b/mex/sources/korderpert/src/fs2000k.mod @@ -1,107 +1,112 @@ -// This file replicates the estimation of the CIA model from -// Frank Schorfheide (2000) "Loss function-based evaluation of DSGE models" -// Journal of Applied Econometrics, 15, 645-670. -// the data are the ones provided on Schorfheide's web site with the programs. -// http://www.econ.upenn.edu/~schorf/programs/dsgesel.ZIP -// You need to have fsdat.m in the same directory as this file. -// This file replicates: -// -the posterior mode as computed by Frank's Gauss programs -// -the parameter mean posterior estimates reported in the paper -// -the model probability (harmonic mean) reported in the paper -// This file was tested with dyn_mat_test_0218.zip -// the smooth shocks are probably stil buggy -// -// The equations are taken from J. Nason and T. Cogley (1994) -// "Testing the implications of long-run neutrality for monetary business -// cycle models" Journal of Applied Econometrics, 9, S37-S70. -// Note that there is an initial minus sign missing in equation (A1), p. S63. -// -// Michel Juillard, February 2004 - -var m m_1 P P_1 c e W R k d n l gy_obs gp_obs Y_obs P_obs y dA P2 c2; -varexo e_a e_m; - -parameters alp bet gam mst rho psi del; - -alp = 0.33; -bet = 0.99; -gam = 0.003; -mst = 1.011; -rho = 0.7; -psi = 0.787; -del = 0.02; - -model (use_dll); -dA = exp(gam+e_a); -log(m) = (1-rho)*log(mst) + rho*log(m_1(-1))+e_m; --P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c2(+1)*P2(+1)*m(+1))=0; -W = l/n; --(psi/(1-psi))*(c*P/(1-n))+l/n = 0; -R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W; -1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0; -c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1); -P*c = m; -m-1+d = l; -e = exp(e_a); -y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a)); -gy_obs = dA*y/y(-1); -gp_obs = (P/P_1(-1))*m_1(-1)/dA; -Y_obs/Y_obs(-1) = gy_obs; -P_obs/P_obs(-1) = gp_obs; -P2 = P(+1); -c2 = c(+1); -m_1 = m; -P_1 = P; -end; - -initval; -k = 6; -m = mst; -P = 2.25; -c = 0.45; -e = 1; -W = 4; -R = 1.02; -d = 0.85; -n = 0.19; -l = 0.86; -y = 0.6; -gy_obs = exp(gam); -gp_obs = exp(-gam); -dA = exp(gam); -end; - -shocks; -var e_a; stderr 0.014; -var e_m; stderr 0.005; -end; - -unit_root_vars P_obs Y_obs; - -steady; - -check; - -estimated_params; -alp, beta_pdf, 0.356, 0.02; -bet, beta_pdf, 0.993, 0.002; -gam, normal_pdf, 0.0085, 0.003; -mst, normal_pdf, 1.0002, 0.007; -rho, beta_pdf, 0.129, 0.223; -psi, beta_pdf, 0.65, 0.05; -del, beta_pdf, 0.01, 0.005; -stderr e_a, inv_gamma_pdf, 0.035449, inf; -stderr e_m, inv_gamma_pdf, 0.008862, inf; -end; - -varobs P_obs Y_obs; - -observation_trends; -P_obs (log(mst)-gam); -Y_obs (gam); -end; - -//options_.useAIM = 1; -estimation(datafile=fsdat,nobs=192,loglinear,mh_replic=2000, - mode_compute=4,mh_nblocks=2,mh_drop=0.45,mh_jscale=0.65); - +// This file replicates the estimation of the CIA model from +// Frank Schorfheide (2000) "Loss function-based evaluation of DSGE models" +// Journal of Applied Econometrics, 15, 645-670. +// the data are the ones provided on Schorfheide's web site with the programs. +// http://www.econ.upenn.edu/~schorf/programs/dsgesel.ZIP +// You need to have fsdat.m in the same directory as this file. +// This file replicates: +// -the posterior mode as computed by Frank's Gauss programs +// -the parameter mean posterior estimates reported in the paper +// -the model probability (harmonic mean) reported in the paper +// This file was tested with dyn_mat_test_0218.zip +// the smooth shocks are probably stil buggy +// +// The equations are taken from J. Nason and T. Cogley (1994) +// "Testing the implications of long-run neutrality for monetary business +// cycle models" Journal of Applied Econometrics, 9, S37-S70. +// Note that there is an initial minus sign missing in equation (A1), p. S63. +// +// Michel Juillard, February 2004 +options_.usePartInfo=0; +//var m P c e W R k d n l gy_obs gp_obs Y_obs P_obs y dA P2 c2; +var m P c e W R k d n l gy_obs gp_obs y dA P2 c2; +varexo e_a e_m; + +parameters alp bet gam mst rho psi del; + +alp = 0.33; +bet = 0.99; +gam = 0.003; +mst = 1.011; +rho = 0.7; +psi = 0.787; +del = 0.02; + +model (use_dll); +dA = exp(gam+e_a); +log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m; +-P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c2(+1)*P2(+1)*m(+1))=0; +W = l/n; +-(psi/(1-psi))*(c*P/(1-n))+l/n = 0; +R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W; +1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0; +c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1); +P*c = m; +m-1+d = l; +e = exp(e_a); +y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a)); +gy_obs = dA*y/y(-1); +gp_obs = (P/P(-1))*m(-1)/dA; +//Y_obs/Y_obs(-1) = gy_obs; +//P_obs/P_obs(-1) = gp_obs; +P2 = P(+1); +c2 = c(+1); +end; + +initval; +m = mst; +P = 2.25; +c = 0.45; +e = 1; +W = 4; +R = 1.02; +k = 6; +d = 0.85; +n = 0.19; +l = 0.86; +gy_obs = exp(gam); +gp_obs = exp(-gam); +// Y_obs = 20000; +// P_obs = 51; +y = 0.6; +dA = exp(gam); + P2=P; + c2=c; +end; + +shocks; +var e_a; stderr 0.014; +var e_m; stderr 0.005; +end; + +//unit_root_vars P_obs Y_obs; + +steady(solve_algo = 2); + +check; + +estimated_params; +alp, beta_pdf, 0.356, 0.02; +bet, beta_pdf, 0.993, 0.002; +gam, normal_pdf, 0.0085, 0.003; +mst, normal_pdf, 1.0002, 0.007; +rho, beta_pdf, 0.129, 0.223; +psi, beta_pdf, 0.65, 0.05; +del, beta_pdf, 0.01, 0.005; +stderr e_a, inv_gamma_pdf, 0.035449, inf; +stderr e_m, inv_gamma_pdf, 0.008862, inf; +end; + +//varobs P_obs Y_obs; +varobs gp_obs gy_obs; + +//observation_trends; +//P_obs (log(mst)-gam); +//Y_obs (gam); +//end; + + +//options_.useAIM = 1; +estimation(datafile=fsdat,nobs=192,loglinear,mh_replic=2000, + mode_compute=4,mh_nblocks=2,mh_drop=0.45,mh_jscale=0.65); + diff --git a/mex/sources/korderpert/src/k_ord_dynare.cpp b/mex/sources/korderpert/src/k_ord_dynare.cpp index 5bdf5f69f..7742fff85 100644 --- a/mex/sources/korderpert/src/k_ord_dynare.cpp +++ b/mex/sources/korderpert/src/k_ord_dynare.cpp @@ -236,11 +236,19 @@ void KordpDynare::calcDerivatives(const Vector& yy, const Vector& xx) Vector& out= *(new Vector(nY)); out.zeros(); + const Vector * llxYYp; // getting around the constantness + if ((nJcols - nExog) > yy.length()){ + llxYYp= (LLxSteady( yy)); + } else { + llxYYp= &yy; + } + const Vector & llxYY=*(llxYYp); + #ifdef DEBUG mexPrintf("k_order_dynaare.cpp: Call eval in calcDerivatives\n"); #endif try { - dynamicDLL.eval( yy, xx, //int nb_row_x, + dynamicDLL.eval( llxYY, xx, //int nb_row_x, params, //int it_, out, g1, NULL); } @@ -364,6 +372,38 @@ void KordpDynare::writeModelInfo(Journal& jr) const } **************/ } +/********************************************************* +* LLxSteady() +* returns ySteady extended with leads and lags suitable for +* passing to _dynamic DLL +*************************************************************/ +Vector * KordpDynare::LLxSteady( const Vector& yS){ + if ((nJcols-nExog) == yS.length()) { + mexPrintf("k_ord_dynare.cpp: Warning in LLxSteady: ySteady already. right size"); + return NULL; + } + // create temporary square 2D matrix size nEndo x nEndo (sparse) + // for the lag, current and lead blocks of the jacobian + Vector * llxSteady = new Vector(nJcols-nExog); + try{ + for (int ll_row=0; ll_row< ll_Incidence->nrows(); ll_row++) + { + // 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]; + } + } + } catch (const TLException& e) { + mexPrintf("Caugth TL exception in LLxSteady: "); + e.print(); + return NULL;// 255; + }catch (...){ + mexPrintf(" Error in LLxSteady - wrong index?"); + } + + return llxSteady; +} /************************************ diff --git a/mex/sources/korderpert/src/k_ord_dynare.h b/mex/sources/korderpert/src/k_ord_dynare.h index 0d82bc627..6c32cb23e 100644 --- a/mex/sources/korderpert/src/k_ord_dynare.h +++ b/mex/sources/korderpert/src/k_ord_dynare.h @@ -212,6 +212,8 @@ public: DynamicModel* clone() const {return new KordpDynare(*this);} void ReorderCols(TwoDMatrix * tdx, const int * 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); diff --git a/mex/sources/korderpert/src/k_order_test_main.cpp b/mex/sources/korderpert/src/k_order_test_main.cpp index 163a3010d..b0d0d4122 100644 --- a/mex/sources/korderpert/src/k_order_test_main.cpp +++ b/mex/sources/korderpert/src/k_order_test_main.cpp @@ -65,28 +65,24 @@ int main(int argc, char* argv[]) 0.0, 0.0250e-3}; npar = 2;//(int)mxGetN(mxFldp); TwoDMatrix * vCov = new TwoDMatrix(npar, npar, (d2Dparams)); - - double dYSparams [31]= { // mxGetData(mxFldp); - 1.0110, 2.2582, 5.8012, 1.0000, 1.0000, 0.5808, 1.0110, 2.2582, - 0.4477, 1.0000, 4.5959, 1.0212, 5.8012, 0.8494, 0.1872, 0.8604, - 1.0030, 1.0080, 1.0000, 1.0000, 0.5808, 1.0030, 2.2582, 0.4477, - 1.0110, 2.2582, 0.4477, 1.0000, 0.1872, 2.2582, 0.4477 + double dYSparams [16]= { // 27 mxGetData(mxFldp); + // 1.0110, 2.2582, 5.8012, 0.5808, + 1.0110, 2.2582, 0.4477, 1.0000 + , 4.5959, 1.0212, 5.8012, 0.8494 + , 0.1872, 0.8604, 1.0030, 1.0080 + , 0.5808, 1.0030, 2.2582, 0.4477 + //, 1.0110, 2.2582, 0.4477, 1.0000, 0.1872, 2.2582, 0.4477 }; - - - // 1.0110, 2.2582, 0.4477, 1.0000, 4.5959, 1.0212, 5.8012, 0.8494, -// 0.1872, 0.8604, 1.0030, 1.0080, 1.0000, 1.0000, 0.5808, 1.0030 -// }; - const int nSteady = 31;//29, 16 (int)mxGetM(mxFldp); + const int nSteady = 16;//27 //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" ); - const int nPred = 4;//6 - nBoth (int)mxGetScalar(mxFldp); + const int nPred = 2;//6 - nBoth (int)mxGetScalar(mxFldp); //mxFldp = mxGetField(dr, 0,"nspred" ); - const int nsPred = 6;//(int)mxGetScalar(mxFldp); + const int nsPred = 4;//(int)mxGetScalar(mxFldp); //mxFldp = mxGetField(dr, 0,"nboth" ); const int nBoth = 2;// (int)mxGetScalar(mxFldp); //mxFldp = mxGetField(dr, 0,"nfwrd" ); @@ -97,7 +93,7 @@ int main(int argc, char* argv[]) //mxFldp = mxGetField(M_, 0,"exo_nbr" ); const int nExog =2;// (int)mxGetScalar(mxFldp); //mxFldp = mxGetField(M_, 0,"endo_nbr" ); - const int nEndo = 18;//16(int)mxGetScalar(mxFldp); + const int nEndo = 16;//16(int)mxGetScalar(mxFldp); //mxFldp = mxGetField(M_, 0,"param_nbr" ); const int nPar = 7;//(int)mxGetScalar(mxFldp); // it_ should be set to M_.maximum_lag @@ -106,37 +102,29 @@ int main(int argc, char* argv[]) int var_order[]//[18] = { - 5, 6, 8, 10, 11, 12, 16, 7, 13, 14, 15, 1, 2, 3, 4, 9, 17, 18 + 5, 6, 8, 10, 11, 12, 14, 7, 13, 1, 2, 3, 4, 9, 15, 16 +// 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 + 1, 5, 21 + , 2, 6, 22 + , 0, 7, 23 + , 0, 8, 24 + , 0, 9, 0 + , 0, 10, 0 + , 3, 11, 0 + , 0, 12, 0 + , 0, 13, 25 + , 0, 14, 0 + , 0, 15, 0 + , 0, 16, 0 + , 4, 17, 0 + , 0, 18, 0 + , 0, 19, 26 + , 0, 20, 27 }; TwoDMatrix * llincidence = new TwoDMatrix ( 3, nEndo, ll_incidence ); @@ -145,9 +133,9 @@ int main(int argc, char* argv[]) 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 nendo = 16;//16(int)mxGetM(mxFldp); const int widthEndo = 6;// (int)mxGetN(mxFldp); - const char * cNamesCharStr= "mPceWRkdnlggYPydPc yp__ A22 __oo oobb bbss ss "; + const char * cNamesCharStr= "mPceWRkdnlggydPc yp A22 __ oo bb ss "; // const char** endoNamesMX= DynareMxArrayToString( mxFldp,nendo,widthEndo); const char** endoNamesMX= DynareMxArrayToString( cNamesCharStr, nendo, widthEndo); #ifdef DEBUG