diff --git a/matlab/dyn_risky_steadystate_solver.m b/matlab/dyn_risky_steadystate_solver.m
index fc559af49bff666afacef8668b7afd3e04e7922b..dd6ad0832298b2842ac5336763e75f99fd37ae85 100644
--- a/matlab/dyn_risky_steadystate_solver.m
+++ b/matlab/dyn_risky_steadystate_solver.m
@@ -86,8 +86,6 @@ order_var = dr.order_var;
endo_nbr = M.endo_nbr;
exo_nbr = M.exo_nbr;
-M.var_order_endo_names = M.endo_names(dr.order_var);
-
[~,dr.i_fwrd_g,i_fwrd_f] = find(lead_lag_incidence(3,order_var));
dr.i_fwrd_f = i_fwrd_f;
nd = nnz(lead_lag_incidence) + M.exo_nbr;
@@ -508,7 +506,6 @@ M_np.endo_names = M.endo_names(v_np);
dr_np = struct();
dr_np = set_state_space(dr_np,M_np,options);
pm.dr_np = dr_np;
-M_np.var_order_endo_names = M_np.endo_names(dr_np.order_var);
pm.M_np = M_np;
pm.i_fwrd_g = find(lead_lag_incidence_np(lead_index,dr_np.order_var)');
diff --git a/matlab/k_order_pert.m b/matlab/k_order_pert.m
index 0fb758448a29117f4166b99b9d06c459bd39bae9..44e40c53ea866caa5b556304311ae76a92736000 100644
--- a/matlab/k_order_pert.m
+++ b/matlab/k_order_pert.m
@@ -20,8 +20,6 @@ function [dr,info] = k_order_pert(dr,M,options)
info = 0;
-M.var_order_endo_names = M.endo_names(dr.order_var);
-
order = options.order;
if order>1 && options.loglinear
diff --git a/mex/sources/k_order_perturbation/k_ord_dynare.cc b/mex/sources/k_order_perturbation/k_ord_dynare.cc
index 744bc0519f2635b2d547da0397b726a42852ee6c..a4f5c5104addb7170f899cd15eb05a6a02a249e0 100644
--- a/mex/sources/k_order_perturbation/k_ord_dynare.cc
+++ b/mex/sources/k_order_perturbation/k_ord_dynare.cc
@@ -25,14 +25,13 @@
#include <utility>
-
KordpDynare::KordpDynare(const std::vector<std::string> &endo,
const std::vector<std::string> &exo, int nexog, int npar,
Vector &ysteady, TwoDMatrix &vcov, Vector &inParams, int nstat,
int npred, int nforw, int nboth, const Vector &nnzd,
int nsteps, int norder,
Journal &jr, std::unique_ptr<DynamicModelAC> dynamicModelFile_arg,
- const std::vector<int> &var_order, const TwoDMatrix &llincidence,
+ const std::vector<int> &dr_order, const TwoDMatrix &llincidence,
std::unique_ptr<TwoDMatrix> g1_arg, std::unique_ptr<TwoDMatrix> g2_arg,
std::unique_ptr<TwoDMatrix> g3_arg) :
nStat{nstat}, nBoth{nboth}, nPred{npred}, nForw{nforw}, nExog{nexog}, nPar{npar},
@@ -44,7 +43,7 @@ KordpDynare::KordpDynare(const std::vector<std::string> &endo,
g1p{std::move(g1_arg)}, g2p{std::move(g2_arg)}, g3p{std::move(g3_arg)},
dynamicModelFile{std::move(dynamicModelFile_arg)}
{
- computeJacobianPermutation(var_order);
+ computeJacobianPermutation(dr_order);
}
void
@@ -195,63 +194,56 @@ KordpDynare::LLxSteady(const Vector &yS, Vector &llxSteady)
llxSteady[static_cast<int>(ll_Incidence.get(ll_row, i))-1] = yS[i];
}
-/************************************
- * Reorder DynareJacobianIndices of variables in a vector according to
- * given int * varOrder together with lead & lag incidence matrix and
- * any the extra columns for exogenous vars, and then,
- * reorders its blocks given by the varOrder and the Dynare++ expectations:
-
- * extra nboth+ npred (t-1) lags
- * varOrder
- static:
- pred
- both
- forward
- * extra both + nforw (t+1) leads, and
- * extra exogen
-
- * so to match the jacobian organisation expected by the Appoximation class
- both + nforw (t+1) leads
- static
- pred
- both
- forward
- nboth+ npred (t-1) lags
- exogen
-************************************/
-
+/*
+ Computes mapping between Dynare and Dynare++ orderings of the (dynamic)
+ variable indices in derivatives.
+
+ If one defines:
+ – y (resp. x) as the vector of all endogenous (size nY), in DR-order (resp.
+ declaration order)
+ – y⁻ (resp. x⁻) as the vector of endogenous that appear at previous period (size nYs),
+ in DR-order (resp. declaration order)
+ – y⁺ (resp. x⁺) as the vector of endogenous that appear at future period (size nYss) in
+ DR-order (resp. declaration order)
+ – u as the vector of exogenous (size nExog)
+
+ In Dynare, the ordering is (x⁻, x, x⁺, u).
+ In Dynare++, the ordering is (y⁺, y, y⁻, u).
+
+ dr_order is typically equal to M_.order_var.
+*/
void
-KordpDynare::computeJacobianPermutation(const std::vector<int> &var_order)
+KordpDynare::computeJacobianPermutation(const std::vector<int> &dr_order)
{
+ // Compute restricted inverse DR-orderings: x⁻→y⁻ and x⁺→y⁺
+ std::vector<int> dr_inv_order_forw(nBoth+nForw), dr_inv_order_pred(nBoth+nPred);
+ std::iota(dr_inv_order_forw.begin(), dr_inv_order_forw.end(), 0);
+ std::sort(dr_inv_order_forw.begin(), dr_inv_order_forw.end(),
+ [&](int i, int j) { return dr_order[nStat+nPred+i] < dr_order[nStat+nPred+j]; });
+ std::iota(dr_inv_order_pred.begin(), dr_inv_order_pred.end(), 0);
+ std::sort(dr_inv_order_pred.begin(), dr_inv_order_pred.end(),
+ [&](int i, int j) { return dr_order[nStat+i] < dr_order[nStat+j]; });
+
+ // Compute restricted DR-orderings: y⁻→x⁻ and y⁺→x⁺
+ std::vector<int> dr_order_forw(nBoth+nForw), dr_order_pred(nBoth+nPred);
+ for (int i = 0; i < nBoth+nForw; i++)
+ dr_order_forw[dr_inv_order_forw[i]] = i;
+ for (int i = 0; i < nBoth+nPred; i++)
+ dr_order_pred[dr_inv_order_pred[i]] = i;
+
+ // Compute Dynare++ → Dynare ordering
dynppToDyn.resize(nJcols);
- // create temporary square 2D matrix size nEndo×nEndo (sparse)
- // for the lag, current and lead blocks of the jacobian
- std::vector<int> tmp(nY);
- int rjoff = nJcols-nExog-1;
-
- for (int ll_row = 0; ll_row < ll_Incidence.nrows(); ll_row++)
- {
- // reorder in orde-var order & populate temporary nEndo (sparse) vector with
- // the lag, current and lead blocks of the jacobian respectively
- for (int i = 0; i < nY; i++)
- tmp[i] = static_cast<int>(ll_Incidence.get(ll_row, var_order[i]-1));
- // write the reordered blocks back to the jacobian
- // in reverse order
- for (int j = nY-1; j >= 0; j--)
- if (tmp[j])
- {
- dynppToDyn[rjoff] = tmp[j]-1;
- rjoff--;
- if (rjoff < 0)
- break;
- }
- }
-
- //add the indices for the nExog exogenous jacobians
- for (int j = nJcols-nExog; j < nJcols; j++)
- dynppToDyn[j] = j;
-
- // Compute reverse ordering
+ int j = 0;
+ for (; j < nYss; j++)
+ dynppToDyn[j] = dr_order_forw[j]+nYs+nY; // Forward variables
+ for (; j < nYss+nY; j++)
+ dynppToDyn[j] = dr_order[j-nYss]+nYs; // Variables in current period
+ for (; j < nYss+nY+nYs; j++)
+ dynppToDyn[j] = dr_order_pred[j-nY-nYss]; // Predetermined variables
+ for (; j < nJcols; j++)
+ dynppToDyn[j] = j; // Exogenous
+
+ // Compute Dynare → Dynare++ ordering
dynToDynpp.resize(nJcols);
for (int i = 0; i < nJcols; i++)
dynToDynpp[dynppToDyn[i]] = i;
diff --git a/mex/sources/k_order_perturbation/k_order_perturbation.cc b/mex/sources/k_order_perturbation/k_order_perturbation.cc
index 438eba56fff4449d4582d7a8bcb494c495547644..5342ba419817ab0a95c4d263af2df7665996d395 100644
--- a/mex/sources/k_order_perturbation/k_order_perturbation.cc
+++ b/mex/sources/k_order_perturbation/k_order_perturbation.cc
@@ -136,14 +136,12 @@ extern "C" {
const int nPar = static_cast<int>(mxGetScalar(mxFldp));
mxFldp = mxGetField(dr, 0, "order_var");
- auto dparams = mxGetPr(mxFldp);
npar = static_cast<int>(mxGetM(mxFldp));
if (npar != nEndo)
DYN_MEX_FUNC_ERR_MSG_TXT("Incorrect number of input var_order vars.");
-
- std::vector<int> var_order_vp(nEndo);
- for (int v = 0; v < nEndo; v++)
- var_order_vp[v] = static_cast<int>(*(dparams++));
+ std::vector<int> dr_order(nEndo);
+ std::transform(mxGetPr(mxFldp), mxGetPr(mxFldp)+npar, dr_order.begin(),
+ [](double x) { return static_cast<int>(x)-1; });
// the lag, current and lead blocks of the jacobian respectively
mxFldp = mxGetField(M_, 0, "lead_lag_incidence");
@@ -160,7 +158,7 @@ extern "C" {
if (NNZD[kOrder-1] == -1)
DYN_MEX_FUNC_ERR_MSG_TXT("The derivatives were not computed for the required order. Make sure that you used the right order option inside the 'stoch_simul' command");
- mxFldp = mxGetField(M_, 0, "var_order_endo_names");
+ mxFldp = mxGetField(M_, 0, "endo_names");
std::vector<std::string> endoNames = DynareMxArrayToString(mxFldp);
mxFldp = mxGetField(M_, 0, "exo_names");
@@ -212,7 +210,7 @@ extern "C" {
KordpDynare dynare(endoNames, exoNames, nExog, nPar,
ySteady, vCov, modParams, nStat, nPred, nForw, nBoth,
NNZD, nSteps, kOrder, journal, std::move(dynamicModelFile),
- var_order_vp, llincidence,
+ dr_order, llincidence,
std::move(g1m), std::move(g2m), std::move(g3m));
// construct main K-order approximation class