Skip to content

GitLab

Commit 79152b55 authored by Johannes Pfeifer's avatar Johannes Pfeifer
Browse files

Merge branch 'master' of ssh://kirikou.dynare.org/srv/d_kirikou/git/dynare

parents a615f025 9d224dba
Hide whitespace changes
Inline Side-by-side
dynare++/tl/cc/t_container.hweb
View file @ 79152b55
......@@ -268,7 +268,7 @@ int getMaxDim() const
@<|TensorContainer::print| code@>=
void print() const
{
printf("Tensor container: nvars=%d, tensors=%d\n", n, m.size());
printf("Tensor container: nvars=%d, tensors=%D\n", n, m.size());
for (const_iterator it = m.begin(); it != m.end(); ++it) {
printf("Symmetry: ");
(*it).first.print();
......
matlab/k_order_pert.m
View file @ 79152b55
......@@ -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
matlab/simult_.m
View file @ 79152b55
......@@ -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
mex/sources/k_order_perturbation/k_order_perturbation.cc
View file @ 79152b55
......@@ -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
preprocessor/ComputingTasks.cc
View file @ 79152b55
......@@ -135,14 +135,6 @@ StochSimulStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsoli
|| mod_file_struct.order_option >= 3)
mod_file_struct.k_order_solver = true;
// Check that option "pruning" is not used with k-order
it = options_list.num_options.find("pruning");
if ((it != options_list.num_options.end() && it->second == "1")
&& mod_file_struct.k_order_solver)
{
cerr << "ERROR: in 'stoch_simul', you cannot use option 'pruning' with 'k_order_solver' option or with 3rd order approximation" << endl;
exit(EXIT_FAILURE);
}
}
void
......
tests/Makefile.am
View file @ 79152b55
......@@ -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 \
......
tests/k_order_perturbation/fs2000k3_p.mod 0 → 100644
View file @ 79152b55
// 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);
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