Commit 3165e2a6 authored by george's avatar george
Browse files

Update, The exit from the walkstochsteady, i.e. the code in the main()...

Update, The exit from the walkstochsteady, i.e. the code in the main() following exit is not complete - it is there only for debugging up to that point.



git-svn-id: https://www.dynare.org/svn/dynare/trunk@2375 ac1d8469-bf42-47a9-8791-bf33cf982152
parent 44a9cd73
// 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);
......@@ -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 <model>_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;i<nY;i++){
if(ll_Incidence->get(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;
}
/************************************
......
......@@ -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);
......
......@@ -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
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment