From 2244af9c627ba2db1d8f906dba7c87933740bbaf Mon Sep 17 00:00:00 2001
From: Michel Juillard <michel.juillard@mjui.fr>
Date: Tue, 31 Jul 2012 21:50:28 +0200
Subject: [PATCH] adding pruning option for 3rd order
---
matlab/k_order_pert.m | 96 ++++++++++++++++++-
matlab/simult_.m | 68 +++++++++++--
.../k_order_perturbation.cc | 65 +++++++++----
tests/Makefile.am | 1 +
tests/k_order_perturbation/fs2000k3_p.mod | 58 +++++++++++
5 files changed, 259 insertions(+), 29 deletions(-)
create mode 100644 tests/k_order_perturbation/fs2000k3_p.mod
diff --git a/matlab/k_order_pert.m b/matlab/k_order_pert.m
index 09eeb955b9..fe0973964e 100644
--- a/matlab/k_order_pert.m
+++ b/matlab/k_order_pert.m
@@ -23,6 +23,9 @@ info = 0;
M.var_order_endo_names = M.endo_names(dr.order_var,:);
order = options.order;
+endo_nbr = M.endo_nbr;
+exo_nbr = M.exo_nbr;
+npred = dr.npred;
switch(order)
case 1
@@ -36,7 +39,25 @@ switch(order)
dr.g_1 = g_1;
dr.g_2 = g_2;
case 3
- [err, g_0, g_1, g_2, g_3] = k_order_perturbation(dr,M,options);
+ if options.pruning
+ [err, g_0, g_1, g_2, g_3, derivs] = k_order_perturbation(dr, ...
+ M,options);
+ dr.ghx = derivs.gy;
+ dr.ghu = derivs.gu;
+ dr.ghxx = unfold2(derivs.gyy,npred);
+ dr.ghxu = derivs.gyu;
+ dr.ghuu = unfold2(derivs.guu,exo_nbr);
+ dr.ghs2 = derivs.gss;
+ dr.ghxxx = unfold3(derivs.gyyy,npred);
+ dr.ghxxu = unfold21(derivs.gyyu,npred,exo_nbr);
+ dr.ghxuu = unfold12(derivs.gyuu,npred,exo_nbr);
+ dr.ghuuu = unfold3(derivs.guuu,exo_nbr);
+ dr.ghxss = derivs.gyss;
+ dr.ghuss = derivs.guss;
+ else
+ [err, g_0, g_1, g_2, g_3] = k_order_perturbation(dr, ...
+ M,options);
+ end
mexErrCheck('k_order_perturbation', err);
dr.g_0 = g_0;
dr.g_1 = g_1;
@@ -46,10 +67,14 @@ switch(order)
error('order > 3 isn''t implemented')
end
+if options.pruning
+ return
+end
+
npred = dr.npred;
-dr.ghx = g_1(:,1:npred);
-dr.ghu = g_1(:,npred+1:end);
+dr.ghx = dr.g_1(:,1:npred);
+dr.ghu = dr.g_1(:,npred+1:end);
if options.loglinear == 1
k = find(dr.kstate(:,2) <= M.maximum_endo_lag+1);
@@ -63,8 +88,6 @@ end
if order > 1
dr.ghs2 = 2*g_0;
- endo_nbr = M.endo_nbr;
- exo_nbr = M.exo_nbr;
s0 = 0;
s1 = 0;
ghxx=zeros(endo_nbr, npred^2);
@@ -97,3 +120,66 @@ if order > 1
dr.ghuu = ghuu;
end
+function y = unfold2(x,n)
+y=zeros(size(x,1),n*n);
+m = 1;
+for i=1:n
+ for j=i:n
+ y(:,(i-1)*n+j)=x(:,m);
+ if j ~= i
+ y(:,(j-1)*n+i)=x(:,m);
+ end
+ m = m+1;
+ end
+end
+
+function y = unfold3(x,n)
+y = zeros(size(x,1),n*n*n);
+m = 1;
+for i=1:n
+ for j=i:n
+ for k=j:n
+ xx = x(:,m);
+ y(:,(i-1)*n*n+(j-1)*n+k) = xx;
+ y(:,(i-1)*n*n+(k-1)*n+j) = xx;
+ y(:,(j-1)*n*n+(k-1)*n+i) = xx;
+ y(:,(i-1)*n*n+(k-1)*n+j) = xx;
+ y(:,(k-1)*n*n+(i-1)*n+j) = xx;
+ y(:,(k-1)*n*n+(j-1)*n+i) = xx;
+ end
+ end
+end
+
+function y = unfold21(x,n1,n2)
+y = zeros(size(x,1),n1*n1*n2);
+m = 1;
+for i=1:n1
+ for j=i:n1
+ for k=1:n2
+ xx = x(:,m);
+ y(:,(i-1)*n1*n2+(j-1)*n2+k) = xx;
+ if j ~= i
+ y(:,(j-1)*n1*n2+(i-1)*n2+i) = xx;
+ end
+ end
+ end
+end
+
+function y = unfold12(x,n1,n2)
+y = zeros(size(x,1),n1*n2*n2);
+m = 1;
+for i=1:n1
+ for j=1:n2
+ for k=j:n2
+ xx = x(:,m);
+ y(:,(i-1)*n2*n2+(j-1)*n2+k) = xx;
+ if k ~= j
+ y(:,(i-1)*n2*n2+(k-1)*n2+j) = xx;
+ end
+ end
+ end
+end
+
+
+
+
\ No newline at end of file
diff --git a/matlab/simult_.m b/matlab/simult_.m
index 1bcc12725b..4e5f23ee90 100644
--- a/matlab/simult_.m
+++ b/matlab/simult_.m
@@ -35,6 +35,8 @@ function y_=simult_(y0,dr,ex_,iorder)
global M_ options_
iter = size(ex_,1);
+endo_nbr = M_.endo_nbr;
+exo_nbr = M_.exo_nbr;
y_ = zeros(size(y0,1),iter+M_.maximum_lag);
y_(:,1) = y0;
@@ -51,19 +53,19 @@ if ~options_.k_order_solver
end
end
-if options_.k_order_solver% Call dynare++ routines.
- ex_ = [zeros(1,M_.exo_nbr); ex_];
+if options_.k_order_solver && ~options_.pruning % Call dynare++ routines.
+ ex_ = [zeros(1,exo_nbr); ex_];
switch options_.order
case 1
- [err, y_] = dynare_simul_(1,dr.nstatic,dr.npred-dr.nboth,dr.nboth,dr.nfwrd,M_.exo_nbr, ...
+ [err, y_] = dynare_simul_(1,dr.nstatic,dr.npred-dr.nboth,dr.nboth,dr.nfwrd,exo_nbr, ...
y_(dr.order_var,1),ex_',M_.Sigma_e,options_.DynareRandomStreams.seed,dr.ys(dr.order_var),...
- zeros(M_.endo_nbr,1),dr.g_1);
+ zeros(endo_nbr,1),dr.g_1);
case 2
- [err, y_] = dynare_simul_(2,dr.nstatic,dr.npred-dr.nboth,dr.nboth,dr.nfwrd,M_.exo_nbr, ...
+ [err, y_] = dynare_simul_(2,dr.nstatic,dr.npred-dr.nboth,dr.nboth,dr.nfwrd,exo_nbr, ...
y_(dr.order_var,1),ex_',M_.Sigma_e,options_.DynareRandomStreams.seed,dr.ys(dr.order_var),dr.g_0, ...
dr.g_1,dr.g_2);
case 3
- [err, y_] = dynare_simul_(3,dr.nstatic,dr.npred-dr.nboth,dr.nboth,dr.nfwrd,M_.exo_nbr, ...
+ [err, y_] = dynare_simul_(3,dr.nstatic,dr.npred-dr.nboth,dr.nboth,dr.nfwrd,exo_nbr, ...
y_(dr.order_var,1),ex_',M_.Sigma_e,options_.DynareRandomStreams.seed,dr.ys(dr.order_var),dr.g_0, ...
dr.g_1,dr.g_2,dr.g_3);
otherwise
@@ -78,11 +80,11 @@ else
else
k2 = [];
end;
- order_var = 1:M_.endo_nbr;
+ order_var = 1:endo_nbr;
dr.order_var = order_var;
else
k2 = dr.kstate(find(dr.kstate(:,2) <= M_.maximum_lag+1),[1 2]);
- k2 = k2(:,1)+(M_.maximum_lag+1-k2(:,2))*M_.endo_nbr;
+ k2 = k2(:,1)+(M_.maximum_lag+1-k2(:,2))*endo_nbr;
order_var = dr.order_var;
end;
@@ -134,6 +136,56 @@ else
y_(dr.order_var,i) = constant + dr.ghx*yhat + dr.ghu*epsilon ...
+ abcOut1 + abcOut2 + abcOut3;
end
+ end
+ case 3
+ % only with pruning
+ ghx = dr.ghx;
+ ghu = dr.ghu;
+ ghxx = dr.ghxx;
+ ghxu = dr.ghxu;
+ ghuu = dr.ghuu;
+ ghs2 = dr.ghs2;
+ ghxxx = dr.ghxxx;
+ ghxxu = dr.ghxxu;
+ ghxuu = dr.ghxuu;
+ ghuuu = dr.ghuuu;
+ ghxss = dr.ghxss;
+ ghuss = dr.ghuss;
+ threads = options_.threads.kronecker.A_times_B_kronecker_C;
+ npred = dr.npred;
+ ipred = dr.nstatic+(1:npred);
+ yhat1 = y0(order_var(k2))-dr.ys(order_var(k2));
+ yhat2 = zeros(npred,1);
+ yhat3 = zeros(npred,1);
+ for i=2:iter+M_.maximum_lag
+ u = ex_(i-1,:)';
+ [gyy, err] = A_times_B_kronecker_C(ghxx,yhat1,threads);
+ mexErrCheck('A_times_B_kronecker_C', err);
+ [guu, err] = A_times_B_kronecker_C(ghuu,u,threads);
+ mexErrCheck('A_times_B_kronecker_C', err);
+ [gyu, err] = A_times_B_kronecker_C(ghxu,yhat1,u,threads);
+ mexErrCheck('A_times_B_kronecker_C', err);
+ y2a = kron(yhat1,yhat1);
+ [gyyy, err] = A_times_B_kronecker_C(ghxxx,y2a,yhat1,threads);
+ mexErrCheck('A_times_B_kronecker_C', err);
+ u2a = kron(u,u);
+ [guuu, err] = A_times_B_kronecker_C(ghuuu,u2a,u,threads);
+ mexErrCheck('A_times_B_kronecker_C', err);
+ yu = kron(yhat1,u);
+ [gyyu, err] = A_times_B_kronecker_C(ghxxu,yhat1,yu,threads);
+ mexErrCheck('A_times_B_kronecker_C', err);
+ [gyuu, err] = A_times_B_kronecker_C(ghxuu,yu,u,threads);
+ mexErrCheck('A_times_B_kronecker_C', err);
+ [gyy12, err] = A_times_B_kronecker_C(ghxx,yhat1,yhat2,threads);
+ mexErrCheck('A_times_B_kronecker_C', err);
+ yhat3 = ghx*yhat3 + gyyy + guuu + 3*(gyyu + gyuu ...
+ + gyy12 + ghxss*yhat1 + ghuss*u);
+ yhat2 = ghx*yhat2 + gyy + guu + 2*gyu + ghs2;
+ yhat1 = ghx*yhat1 + ghu*u;
+ y_(order_var,i) = yhat1 + (1/2)*yhat2 + (1/6)*yhat3;
+ yhat1 = yhat1(ipred);
+ yhat2 = yhat2(ipred);
+ yhat3 = yhat3(ipred);
end
end
end
\ No newline at end of file
diff --git a/mex/sources/k_order_perturbation/k_order_perturbation.cc b/mex/sources/k_order_perturbation/k_order_perturbation.cc
index fef539f19e..59100d23c9 100644
--- a/mex/sources/k_order_perturbation/k_order_perturbation.cc
+++ b/mex/sources/k_order_perturbation/k_order_perturbation.cc
@@ -62,6 +62,16 @@ DynareMxArrayToString(const mxArray *mxFldp, const int len, const int width, vec
out[j] += cNamesCharStr[j+i*len];
}
+void copy_derivatives(mxArray *destin,const Symmetry &sym,const FGSContainer *derivs,const std::string &fieldname)
+{
+ const TwoDMatrix* x = derivs->get(sym);
+ int n = x->numRows();
+ int m = x->numCols();
+ mxArray *tmp = mxCreateDoubleMatrix(n, m, mxREAL);
+ memcpy(mxGetPr(tmp),x->getData().base(),n*m*sizeof(double));
+ mxSetField(destin,0,fieldname.c_str(),tmp);
+}
+
extern "C" {
void
@@ -89,10 +99,10 @@ extern "C" {
else
kOrder = 1;
- if (kOrder == 1 && nlhs != 2)
- DYN_MEX_FUNC_ERR_MSG_TXT("k_order_perturbation at order 1 requires exactly 2 arguments in output");
- else if (kOrder > 1 && nlhs != kOrder+2)
- DYN_MEX_FUNC_ERR_MSG_TXT("k_order_perturbation at order > 1 requires exactly order+2 arguments in output");
+ //if (kOrder == 1 && nlhs != 2)
+ // DYN_MEX_FUNC_ERR_MSG_TXT("k_order_perturbation at order 1 requires exactly 2 arguments in output");
+ //else if (kOrder > 1 && nlhs != kOrder+2)
+ // DYN_MEX_FUNC_ERR_MSG_TXT("k_order_perturbation at order > 1 requires exactly order+2 arguments in output");
double qz_criterium = 1+1e-6;
mxFldp = mxGetField(options_, 0, "qz_criterium");
@@ -250,6 +260,8 @@ extern "C" {
// run stochastic steady
app.walkStochSteady();
+ app.get_rule_ders()->print();
+
/* Write derivative outputs into memory map */
map<string, ConstTwoDMatrix> mm;
app.getFoldDecisionRule().writeMMap(mm, string());
@@ -275,18 +287,40 @@ extern "C" {
for (map<string, ConstTwoDMatrix>::const_iterator cit = mm.begin();
((cit != mm.end()) && (ii < nlhs)); ++cit)
{
- {
- plhs[ii] = mxCreateDoubleMatrix((*cit).second.numRows(), (*cit).second.numCols(), mxREAL);
-
- // Copy Dynare++ matrix into MATLAB matrix
- const ConstVector &vec = (*cit).second.getData();
- assert(vec.skip() == 1);
- memcpy(mxGetPr(plhs[ii]), vec.base(), vec.length() * sizeof(double));
-
- ++ii;
- }
+ plhs[ii] = mxCreateDoubleMatrix((*cit).second.numRows(), (*cit).second.numCols(), mxREAL);
+
+ // Copy Dynare++ matrix into MATLAB matrix
+ const ConstVector &vec = (*cit).second.getData();
+ assert(vec.skip() == 1);
+ memcpy(mxGetPr(plhs[ii]), vec.base(), vec.length() * sizeof(double));
+
+ ++ii;
+
}
- }
+ if (kOrder == 3 && nlhs > 4)
+ {
+ const FGSContainer *derivs = app.get_rule_ders();
+ const std::string fieldnames[] = {"gy", "gu", "gyy", "gyu", "guu", "gss",
+ "gyyy", "gyyu", "gyuu", "guuu", "gyss", "guss"};
+ // creates the char** expected by mxCreateStructMatrix()
+ const char* c_fieldnames[12];
+ for (int i=0; i < 12;++i)
+ c_fieldnames[i] = fieldnames[i].c_str();
+ plhs[ii] = mxCreateStructMatrix(1,1,12,c_fieldnames);
+ copy_derivatives(plhs[ii],Symmetry(1,0,0,0),derivs,"gy");
+ copy_derivatives(plhs[ii],Symmetry(0,1,0,0),derivs,"gu");
+ copy_derivatives(plhs[ii],Symmetry(2,0,0,0),derivs,"gyy");
+ copy_derivatives(plhs[ii],Symmetry(0,2,0,0),derivs,"guu");
+ copy_derivatives(plhs[ii],Symmetry(1,1,0,0),derivs,"gyu");
+ copy_derivatives(plhs[ii],Symmetry(0,0,0,2),derivs,"gss");
+ copy_derivatives(plhs[ii],Symmetry(3,0,0,0),derivs,"gyyy");
+ copy_derivatives(plhs[ii],Symmetry(0,3,0,0),derivs,"guuu");
+ copy_derivatives(plhs[ii],Symmetry(2,1,0,0),derivs,"gyyu");
+ copy_derivatives(plhs[ii],Symmetry(1,2,0,0),derivs,"gyuu");
+ copy_derivatives(plhs[ii],Symmetry(1,0,0,2),derivs,"gyss");
+ copy_derivatives(plhs[ii],Symmetry(0,1,0,2),derivs,"guss");
+ }
+ }
}
catch (const KordException &e)
{
@@ -320,5 +354,4 @@ extern "C" {
plhs[0] = mxCreateDoubleScalar(0);
} // end of mexFunction()
} // end of extern C
-
#endif // ifdef MATLAB_MEX_FILE to exclude mexFunction for other applications
diff --git a/tests/Makefile.am b/tests/Makefile.am
index e76a91a592..e1f8792521 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -49,6 +49,7 @@ MODFILES = \
k_order_perturbation/fs2000k2_m.mod \
k_order_perturbation/fs2000k_1_m.mod \
k_order_perturbation/fs2000k3_m.mod \
+ k_order_perturbation/fs2000k3_p.mod \
partial_information/PItest3aHc0PCLsimModPiYrVarobsAll.mod \
partial_information/PItest3aHc0PCLsimModPiYrVarobsCNR.mod \
arima/mod1.mod \
diff --git a/tests/k_order_perturbation/fs2000k3_p.mod b/tests/k_order_perturbation/fs2000k3_p.mod
new file mode 100644
index 0000000000..2d0c6bff12
--- /dev/null
+++ b/tests/k_order_perturbation/fs2000k3_p.mod
@@ -0,0 +1,58 @@
+// checks whether second order coefficients are the same with order=2 and order=3 with k_order_solver=1
+
+var m P c e W R k d n l gy_obs gp_obs y dA ;
+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;
+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)))))/(c(+2)*P(+2)*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;
+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;
+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;
+
+steady;
+
+stoch_simul(order=3,periods=1000,drop=0,pruning);
+
--
GitLab