diff --git a/matlab/dyn_risky_steadystate_solver.m b/matlab/dyn_risky_steadystate_solver.m
index 04c50b84b234287a994130c8e792906f8119f810..fc559af49bff666afacef8668b7afd3e04e7922b 100644
--- a/matlab/dyn_risky_steadystate_solver.m
+++ b/matlab/dyn_risky_steadystate_solver.m
@@ -346,8 +346,12 @@ if nargout > 1
nu2 = exo_nbr*(exo_nbr+1)/2;
nu3 = exo_nbr*(exo_nbr+1)*(exo_nbr+2)/3;
M_np.NZZDerivatives = [nnz(d1_np); nnz(d2_np); nnz(d3_np)];
- [err,g_0, g_1, g_2, g_3] = k_order_perturbation(dr_np,M_np,options,d1_np,d2_np,d3_np);
+ [err, dynpp_derivs] = k_order_perturbation(dr_np,M_np,options,d1_np,d2_np,d3_np);
mexErrCheck('k_order_perturbation', err);
+ g_0 = dynpp_derivs.g_0;
+ g_1 = dynpp_derivs.g_1;
+ g_2 = dynpp_derivs.g_2;
+ g_3 = dynpp_derivs.g_3;
gu1 = g_1(i_fwrd_g,end-exo_nbr+1:end);
ghuu = unfold2(g_2(:,end-nu2+1:end),exo_nbr);
@@ -547,4 +551,4 @@ for i=1:n
m2 = m2 - (n-i+1);
kk = kk + nx*nx;
end
-end
\ No newline at end of file
+end
diff --git a/matlab/k_order_pert.m b/matlab/k_order_pert.m
index 892cec68b1ae17f12102e6779a0301133b0b6705..3a903b5345f740ab5ac01414839f85548e828647 100644
--- a/matlab/k_order_pert.m
+++ b/matlab/k_order_pert.m
@@ -28,9 +28,6 @@ M.var_order_endo_names = char(M.var_order_endo_names);
M.exo_names = char(M.exo_names);
order = options.order;
-endo_nbr = M.endo_nbr;
-exo_nbr = M.exo_nbr;
-nspred = M.nspred;
if order>1 && options.loglinear
error('The loglinear-option currently only works at order 1')
@@ -40,106 +37,49 @@ if M.maximum_endo_lead == 0 && order>1
'backward models'])
end
-switch(order)
- case 1
- [err, g_1] = k_order_perturbation(dr,M,options);
- if err
- info(1)=9;
- return
- end
- dr.g_1 = g_1;
- case 2
- [err, g_0, g_1, g_2] = k_order_perturbation(dr,M,options);
- if err
- info(1)=9;
- return
- end
- dr.g_0 = g_0;
- dr.g_1 = g_1;
- dr.g_2 = g_2;
- case 3
- if options.pruning
- [err, g_0, g_1, g_2, g_3, derivs] = k_order_perturbation(dr,M,options);
- if err
- info(1)=9;
- return
- end
- else
- [err, g_0, g_1, g_2, g_3] = k_order_perturbation(dr,M,options);
- if err
- info(1)=9;
- return
- end
- end
- dr.g_0 = g_0;
- dr.g_1 = g_1;
- dr.g_2 = g_2;
- dr.g_3 = g_3;
- otherwise
+if order > 3
error('order > 3 isn''t implemented')
end
-% Now fill in dr.ghx, dr.ghu...
+[err, dynpp_derivs, dyn_derivs] = k_order_perturbation(dr,M,options);
+if err
+ info(1)=9;
+ return
+end
+
+% Copy Dynare++ tensors
+
+for i = 0:order
+ gname = [ 'g_' num2str(i) ];
+ dr.(gname) = dynpp_derivs.(gname);
+end
-if options.pruning && order == 3
- dr.ghx = derivs.gy;
- dr.ghu = derivs.gu;
- dr.ghxx = derivs.gyy;
- dr.ghxu = derivs.gyu;
- dr.ghuu = derivs.guu;
- dr.ghs2 = derivs.gss;
- dr.ghxxx = derivs.gyyy;
- dr.ghxxu = derivs.gyyu;
- dr.ghxuu = derivs.gyuu;
- dr.ghuuu = derivs.guuu;
- dr.ghxss = derivs.gyss;
- dr.ghuss = derivs.guss;
-else
- nspred = M.nspred;
+% Fill equivalent Dynare matrices (up to 3rd order)
- dr.ghx = dr.g_1(:,1:nspred);
- dr.ghu = dr.g_1(:,nspred+1:end);
+dr.ghx = dyn_derivs.gy;
+dr.ghu = dyn_derivs.gu;
- if options.loglinear
- k = find(dr.kstate(:,2) <= M.maximum_endo_lag+1);
- klag = dr.kstate(k,[1 2]);
- k1 = dr.order_var;
- dr.ghx = repmat(1./dr.ys(k1),1,size(dr.ghx,2)).*dr.ghx.* ...
- repmat(dr.ys(k1(klag(:,1)))',size(dr.ghx,1),1);
- dr.ghu = repmat(1./dr.ys(k1),1,size(dr.ghu,2)).*dr.ghu;
- end
+if options.loglinear
+ k = find(dr.kstate(:,2) <= M.maximum_endo_lag+1);
+ klag = dr.kstate(k,[1 2]);
+ k1 = dr.order_var;
+ dr.ghx = repmat(1./dr.ys(k1),1,size(dr.ghx,2)).*dr.ghx.* ...
+ repmat(dr.ys(k1(klag(:,1)))',size(dr.ghx,1),1);
+ dr.ghu = repmat(1./dr.ys(k1),1,size(dr.ghu,2)).*dr.ghu;
+end
+
+if order >= 2
+ dr.ghxx = dyn_derivs.gyy;
+ dr.ghxu = dyn_derivs.gyu;
+ dr.ghuu = dyn_derivs.guu;
+ dr.ghs2 = dyn_derivs.gss;
+end
- if order > 1
- dr.ghs2 = 2*g_0;
- s0 = 0;
- s1 = 0;
- ghxx=zeros(endo_nbr, nspred^2);
- ghxu=zeros(endo_nbr, nspred*exo_nbr);
- ghuu=zeros(endo_nbr, exo_nbr^2);
- for i=1:size(g_2,2)
- if s0 < nspred && s1 < nspred
- ghxx(:,s0*nspred+s1+1) = 2*g_2(:,i);
- if s1 > s0
- ghxx(:,s1*nspred+s0+1) = 2*g_2(:,i);
- end
- elseif s0 < nspred && s1 < nspred+exo_nbr
- ghxu(:,(s0*exo_nbr+s1-nspred+1)) = 2*g_2(:,i);
- elseif s0 < nspred+exo_nbr && s1 < nspred+exo_nbr
- ghuu(:,(s0-nspred)*exo_nbr+s1-nspred +1) = 2*g_2(:,i);
- if s1 > s0
- ghuu(:,(s1-nspred)*exo_nbr+s0-nspred+1) = 2*g_2(:,i);
- end
- else
- error('dr1:k_order_perturbation:g_2','Unaccounted columns in g_2');
- end
- s1 = s1+1;
- if s1 == nspred+exo_nbr
- s0 = s0+1;
- s1 = s0;
- end
- end % for loop
- dr.ghxx = ghxx;
- dr.ghxu = ghxu;
- dr.ghuu = ghuu;
- end
+if order >= 3
+ dr.ghxxx = dyn_derivs.gyyy;
+ dr.ghxxu = dyn_derivs.gyyu;
+ dr.ghxuu = dyn_derivs.gyuu;
+ dr.ghuuu = dyn_derivs.guuu;
+ dr.ghxss = dyn_derivs.gyss;
+ dr.ghuss = dyn_derivs.guss;
end
diff --git a/matlab/mex/k_order_perturbation.m b/matlab/mex/k_order_perturbation.m
index 7c10d0efdf86ac44d765b331cf6fc9a86d81f4b7..2f4b1eec96723d045e5b688d6571415b91a2bc81 100644
--- a/matlab/mex/k_order_perturbation.m
+++ b/matlab/mex/k_order_perturbation.m
@@ -1,54 +1,31 @@
-% [err, g_0, g_1, g_2, g_3, derivs] = k_order_perturbation(dr,DynareModel,DynareOptions)
+% [err, dynpp_derivs, dyn_derivs] = k_order_perturbation(dr,DynareModel,DynareOptions,g1,g2,g3)
% computes a k_order_petrubation solution for k=1,2,3
%
% INPUTS
% dr: struct describing the reduced form solution of the model.
% DynareModel: struct jobs's parameters
% DynareOptions: struct job's options
+% g1: matrix First order derivatives of the model (optional)
+% g2: matrix Second order derivatives of the model (optional)
+% g3: matrix Third order derivatives of the model (optional)
%
% OUTPUTS
-% err: double err code (currently unused)
-% g_0: vector dynare++ output. Constant effect of future volatility (in
-% dr.order_var order) on the decision rule
-% (=ghs2/2). Contains zero when order of
-% approximation is 1.
-% g_1: matrix dynare++ output. First order Taylor coefficients
-% of decision rule. When order of approximation
-% is 3, the. Contains both the effect of
-% state endogenous variable and shocks. The rows
-% are in dr.order_var order. The columns are in
-% dr.order_var order of state endogenous
-% variables and shocks
-% g_2: matrix dynare++ output. Second order Taylor coefficients of decision
-% rule. Contains both the effect of state endogenous
-% variable and shocks. The rows are in dr.order_var
-% order. Each row corresponds to the vectorized
-% version of the lower triangle of the Hessian
-% matrix. The Taylor coefficient (1/2) is
-% included. The columns of the Hessian matrix are in
-% dr.order_var order of state endogenous variables
-% and shocks
-% g_3: matrix dynare++ output. Third order Taylor coefficients of decision
-% rule. Contains both the effect of state endogenous
-% variable and shocks. The rows are in dr.order_var
-% order. Each row corresponds to the vectorized
-% version of the 3rd order derivatives tensor where each
-% combination of variables appears only once.
-% The Taylor coefficient (1/6) is
-% included. Inside the tensor, the variables are in
-% dr.order_var order of state endogenous variables
-% and shocks
-% derivs struct contains the original derivatives of the
-% decision function (ghx, ghu, ghxx, ghxu, ghuu,
-% ghs2, ghxxx, ghxxu, ghxuu,ghuuu, ghxss, ghuss),
-% keeping the effect of future volatility
-% separate (in ghs2, ghxss and ghuss). The
-% derivatives matrices contain full versions of
-% the Hessian matrices and 3rd order
-% tensor. Symmetric derivatives are repeated. The
-% Taylor coefficients (1/2 and 1/6) aren't
+% err: double error code
+% dynpp_derivs struct Derivatives of the decision rule in Dynare++ format.
+% The tensors are folded and the Taylor coefficients (1/n!) are
% included.
-% k_order_peturbation is a compiled MEX function. It's source code is in
+% Fieldnames are g_0, g_1, g_2, …
+% dyn_derivs struct Derivatives of the decision rule in Dynare format, up to third order.
+% Matrices are dense and unfolded. The Taylor coefficients (1/2
+% and 1/6) aren't included.
+% The derivatives w.r.t. different categories of variables
+% are separated.
+% Fieldnames are:
+% + gy, gu
+% + if order ≥ 2: gyy, gyu, guu, gss
+% + if order ≥ 3: gyyy, gyyu, gyuu, guuu, gyss, guss
+%
+% k_order_perturbation is a compiled MEX function. It's source code is in
% dynare/mex/sources/k_order_perturbation.cc and it uses code provided by
% dynare++
diff --git a/mex/sources/k_order_perturbation/k_order_perturbation.cc b/mex/sources/k_order_perturbation/k_order_perturbation.cc
index 27da80e9ab6e9f1b59a419a6e5fc6c29e5a8ff70..f518895eba0c0fdfeafdefd90dbb94a77ede38cc 100644
--- a/mex/sources/k_order_perturbation/k_order_perturbation.cc
+++ b/mex/sources/k_order_perturbation/k_order_perturbation.cc
@@ -39,6 +39,13 @@
#include "dynmex.h"
+/* Vector for storing field names like “g_0”, “g_1”, …
+ A static structure is needed since MATLAB apparently does not create its own
+ copy of the strings (contrary to what is said at:
+ https://fr.mathworks.com/matlabcentral/answers/315937-mxcreatestructarray-and-mxcreatestructmatrix-field-name-memory-management
+ ) */
+std::vector<std::string> g_fieldnames;
+
/* Convert MATLAB Dynare endo and exo names array to a vector<string> array of
string pointers. MATLAB “mx” function returns a long string concatenated by
columns rather than rows hence a rather low level approach is needed. */
@@ -57,7 +64,7 @@ DynareMxArrayToString(const mxArray *mxFldp, int len, int width, std::vector<std
}
void
-copy_derivatives(mxArray *destin, const Symmetry &sym, const FGSContainer &derivs, const std::string &fieldname)
+copy_derivatives(mxArray *destin, const Symmetry &sym, const FGSContainer &derivs, const char *fieldname)
{
const FGSTensor &x = derivs.get(sym);
auto x_unfolded = x.unfold();
@@ -65,7 +72,7 @@ copy_derivatives(mxArray *destin, const Symmetry &sym, const FGSContainer &deriv
int m = x_unfolded->numCols();
mxArray *tmp = mxCreateDoubleMatrix(n, m, mxREAL);
std::copy_n(x_unfolded->getData().base(), n*m, mxGetPr(tmp));
- mxSetField(destin, 0, fieldname.c_str(), tmp);
+ mxSetField(destin, 0, fieldname, tmp);
}
extern "C" {
@@ -224,64 +231,57 @@ extern "C" {
// run stochastic steady
app.walkStochSteady();
- /* Write derivative outputs into memory map */
- std::map<std::string, ConstTwoDMatrix> mm;
- app.getFoldDecisionRule().writeMMap(mm, "");
+ const FoldDecisionRule &fdr = app.getFoldDecisionRule();
- // get latest ysteady
- ySteady = dynare.getSteady();
+ // Add possibly missing field names
+ for (int i = static_cast<int>(g_fieldnames.size()); i <= kOrder; i++)
+ g_fieldnames.emplace_back("g_" + std::to_string(i));
+ // Create structure for storing derivatives in Dynare++ format
+ const char *g_fieldnames_c[kOrder+1];
+ for (int i = 0; i <= kOrder; i++)
+ g_fieldnames_c[i] = g_fieldnames[i].c_str();
+ plhs[1] = mxCreateStructMatrix(1, 1, kOrder+1, g_fieldnames_c);
- if (kOrder == 1)
+ // Fill that structure
+ for (int i = 0; i <= kOrder; i++)
{
- /* Set the output pointer to the output matrix ysteady. */
- auto cit = mm.begin();
- ++cit;
- plhs[1] = mxCreateDoubleMatrix(cit->second.numRows(), cit->second.numCols(), mxREAL);
-
- // Copy Dynare++ matrix into MATLAB matrix
- const ConstVector &vec = cit->second.getData();
+ const FFSTensor &t = fdr.get(Symmetry{i});
+ mxArray *tmp = mxCreateDoubleMatrix(t.numRows(), t.numCols(), mxREAL);
+ const ConstVector &vec = t.getData();
assert(vec.skip() == 1);
- std::copy_n(vec.base(), vec.length(), mxGetPr(plhs[1]));
+ std::copy_n(vec.base(), vec.length(), mxGetPr(tmp));
+ mxSetField(plhs[1], 0, ("g_" + std::to_string(i)).c_str(), tmp);
}
- if (kOrder >= 2)
+
+ if (nlhs > 2)
{
- int ii = 1;
- for (auto cit = mm.begin(); cit != mm.end() && ii < nlhs; ++cit)
+ /* Return as 3rd argument a struct containing derivatives in Dynare
+ format (unfolded matrices, without Taylor coefficient) up to 3rd
+ order */
+ const FGSContainer &derivs = app.get_rule_ders();
+
+ size_t nfields = (kOrder == 1 ? 2 : (kOrder == 2 ? 6 : 12));
+ const char *c_fieldnames[] = { "gy", "gu", "gyy", "gyu", "guu", "gss",
+ "gyyy", "gyyu", "gyuu", "guuu", "gyss", "guss" };
+ plhs[2] = mxCreateStructMatrix(1, 1, nfields, c_fieldnames);
+
+ copy_derivatives(plhs[2], Symmetry{1, 0, 0, 0}, derivs, "gy");
+ copy_derivatives(plhs[2], Symmetry{0, 1, 0, 0}, derivs, "gu");
+ if (kOrder >= 2)
{
- 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);
- std::copy_n(vec.base(), vec.length(), mxGetPr(plhs[ii]));
-
- ++ii;
-
+ copy_derivatives(plhs[2], Symmetry{2, 0, 0, 0}, derivs, "gyy");
+ copy_derivatives(plhs[2], Symmetry{0, 2, 0, 0}, derivs, "guu");
+ copy_derivatives(plhs[2], Symmetry{1, 1, 0, 0}, derivs, "gyu");
+ copy_derivatives(plhs[2], Symmetry{0, 0, 0, 2}, derivs, "gss");
}
- if (kOrder == 3 && nlhs > 5)
+ if (kOrder >= 3)
{
- /* Return as 5th argument a struct containing *unfolded* matrices
- for 3rd order decision rule */
- 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");
+ copy_derivatives(plhs[2], Symmetry{3, 0, 0, 0}, derivs, "gyyy");
+ copy_derivatives(plhs[2], Symmetry{0, 3, 0, 0}, derivs, "guuu");
+ copy_derivatives(plhs[2], Symmetry{2, 1, 0, 0}, derivs, "gyyu");
+ copy_derivatives(plhs[2], Symmetry{1, 2, 0, 0}, derivs, "gyuu");
+ copy_derivatives(plhs[2], Symmetry{1, 0, 0, 2}, derivs, "gyss");
+ copy_derivatives(plhs[2], Symmetry{0, 1, 0, 2}, derivs, "guss");
}
}
}