diff --git a/mex/sources/k_order_welfare/approximation_welfare.cc b/mex/sources/k_order_welfare/approximation_welfare.cc index 5e3b527f9..ccd82854c 100644 --- a/mex/sources/k_order_welfare/approximation_welfare.cc +++ b/mex/sources/k_order_welfare/approximation_welfare.cc @@ -24,8 +24,7 @@ #include "approximation_welfare.hh" ApproximationWelfare::ApproximationWelfare(KordwDynare &w, double discount_factor_arg, const FGSContainer &rule_ders_arg, const FGSContainer &rule_ders_s_arg, Journal &j) - : welfare{w}, discount_factor(discount_factor_arg), cond(1), - nvs{welfare.getModel().nys(), welfare.getModel().nexog(), welfare.getModel().nexog(), 1}, journal{j} + : welfare{w}, discount_factor(discount_factor_arg), nvs{welfare.getModel().nys(), welfare.getModel().nexog(), welfare.getModel().nexog(), 1}, journal{j} { rule_ders = std::make_unique(rule_ders_arg); rule_ders_s = std::make_unique(rule_ders_s_arg); @@ -42,23 +41,12 @@ ApproximationWelfare::approxAtSteady() KOrderWelfare korderwel(welfare.getModel().nstat(), welfare.getModel().npred(), welfare.getModel().nboth(), welfare.getModel().nforw(), welfare.getModel().nexog(), welfare.getModel().order(), discount_factor, welfare.getPlannerObjDerivatives(), get_rule_ders(), get_rule_ders_s(), welfare.getModel().getVcov(), journal); for (int k = 1; k <= welfare.getModel().order(); k++) korderwel.performStep(k); - saveRuleDerivs(korderwel.getFoldU(), korderwel.getFoldW()); - // construct the welfare approximations - calcStochShift(); //conditional welfare + saveRuleDerivs(korderwel.getFoldW()); - // construct the resulting decision rules - unc_fdr = std::make_unique(*unc_ders, welfare.getModel().nys(), welfare.getModel().nexog(), welfare.getModel().getSteady()); + // construct the resulting decision rule cond_fdr = std::make_unique(*cond_ders, welfare.getModel().nys(), welfare.getModel().nexog(), welfare.getModel().getSteady()); } -/* This just returns ‘fdr’ with a check that it is created. */ -const FoldDecisionRule & -ApproximationWelfare::getFoldUncWel() const -{ - KORD_RAISE_IF(!unc_fdr, - "Folded decision rule has not been created in ApproximationWelfare::getFoldUncWel"); - return *unc_fdr; -} const FoldDecisionRule & ApproximationWelfare::getFoldCondWel() const { @@ -68,29 +56,7 @@ ApproximationWelfare::getFoldCondWel() const } void -ApproximationWelfare::saveRuleDerivs(const FGSContainer &U, const FGSContainer &W) +ApproximationWelfare::saveRuleDerivs(const FGSContainer &W) { - unc_ders = std::make_unique(U); cond_ders = std::make_unique(W); -} - -/* This method calculates the shift of the conditional welfare function w.r.t its steady-state value because of the presence of stochastic shocks, i.e. - 1 -W ≈ Wbar + ∑ ── W_σᵈ - d=1 d! -*/ -void -ApproximationWelfare::calcStochShift() -{ - cond.zeros(); - cond.add(1.0/(1.0-discount_factor), welfare.getResid()); - KORD_RAISE_IF(cond.length() != 1, - "Wrong length of output vector for ApproximationWelfare::calcStochShift"); - int dfac = 1; - for (int d = 1; d <= cond_ders->getMaxDim(); d++, dfac *= d) - if (KOrderWelfare::is_even(d)) - { - Symmetry sym{0, 0, 0, d}; - cond.add(1.0/dfac, cond_ders->get(sym).getData()); - } -} +} \ No newline at end of file diff --git a/mex/sources/k_order_welfare/approximation_welfare.hh b/mex/sources/k_order_welfare/approximation_welfare.hh index c2d2d166f..9c22eae25 100644 --- a/mex/sources/k_order_welfare/approximation_welfare.hh +++ b/mex/sources/k_order_welfare/approximation_welfare.hh @@ -31,14 +31,9 @@ class ApproximationWelfare double discount_factor; std::unique_ptr rule_ders; std::unique_ptr rule_ders_s; - std::unique_ptr unc_ders; std::unique_ptr cond_ders; - std::unique_ptr unc_fdr; std::unique_ptr cond_fdr; - Vector cond; - // const FNormalMoments mom; IntSequence nvs; - // TwoDMatrix ss; Journal &journal; public: ApproximationWelfare(KordwDynare &w, double discount_factor, const FGSContainer &rule_ders, const FGSContainer &rule_ders_s, Journal &j); @@ -59,30 +54,16 @@ public: return *rule_ders_s; } const FGSContainer & - get_unc_ders() const - { - return *unc_ders; - } - const FGSContainer & get_cond_ders() const { return *cond_ders; } - const Vector & - getCond() const - { - return cond; - } - const FoldDecisionRule & getFoldUncWel() const; - const FoldDecisionRule & getUnfoldUncWel() const; const FoldDecisionRule & getFoldCondWel() const; - const FoldDecisionRule & getUnfoldCondWel() const; void approxAtSteady(); protected: - void saveRuleDerivs(const FGSContainer &U, const FGSContainer &W); - void calcStochShift(); + void saveRuleDerivs(const FGSContainer &W); }; #endif diff --git a/mex/sources/k_order_welfare/k_order_welfare.cc b/mex/sources/k_order_welfare/k_order_welfare.cc index 266d13338..53501411e 100644 --- a/mex/sources/k_order_welfare/k_order_welfare.cc +++ b/mex/sources/k_order_welfare/k_order_welfare.cc @@ -59,7 +59,6 @@ DynareMxArrayToString(const mxArray *mxFldp) ) */ -std::vector U_fieldnames; std::vector W_fieldnames; void @@ -278,113 +277,25 @@ extern "C" { ApproximationWelfare appwel(welfare, discount_factor, app.get_rule_ders(), app.get_rule_ders_s(), journal); appwel.approxAtSteady(); - // Conditional welfare approximation at non-stochastic steady state when stochastic shocks are enabled - const Vector &cond_approx = appwel.getCond(); - plhs[0] = mxCreateDoubleMatrix(cond_approx.length(), 1, mxREAL); - std::copy_n(cond_approx.base(), cond_approx.length(), mxGetPr(plhs[0])); + const FoldDecisionRule &cond_fdr = appwel.getFoldCondWel(); + // Add possibly missing field names + for (int i = static_cast(W_fieldnames.size()); i <= kOrder; i++) + W_fieldnames.emplace_back("W_" + std::to_string(i)); + // Create structure for storing derivatives in Dynare++ format + const char *W_fieldnames_c[kOrder+1]; + for (int i = 0; i <= kOrder; i++) + W_fieldnames_c[i] = W_fieldnames[i].c_str(); + plhs[0] = mxCreateStructMatrix(1, 1, kOrder+1, W_fieldnames_c); - if (nlhs > 1) - { - const FoldDecisionRule &unc_fdr = appwel.getFoldUncWel(); - // Add possibly missing field names - for (int i = static_cast(U_fieldnames.size()); i <= kOrder; i++) - U_fieldnames.emplace_back("U_" + std::to_string(i)); - // Create structure for storing derivatives in Dynare++ format - const char *U_fieldnames_c[kOrder+1]; - for (int i = 0; i <= kOrder; i++) - U_fieldnames_c[i] = U_fieldnames[i].c_str(); - plhs[1] = mxCreateStructMatrix(1, 1, kOrder+1, U_fieldnames_c); - - // Fill that structure - for (int i = 0; i <= kOrder; i++) - { - const FFSTensor &t = unc_fdr.get(Symmetry{i}); - mxArray *tmp = mxCreateDoubleMatrix(t.nrows(), t.ncols(), mxREAL); - const ConstVector &vec = t.getData(); - assert(vec.skip() == 1); - std::copy_n(vec.base(), vec.length(), mxGetPr(tmp)); - mxSetField(plhs[1], 0, ("U_" + std::to_string(i)).c_str(), tmp); - } - - if (nlhs > 2) - { - const FoldDecisionRule &cond_fdr = appwel.getFoldCondWel(); - // Add possibly missing field names - for (int i = static_cast(W_fieldnames.size()); i <= kOrder; i++) - W_fieldnames.emplace_back("W_" + std::to_string(i)); - // Create structure for storing derivatives in Dynare++ format - const char *W_fieldnames_c[kOrder+1]; - for (int i = 0; i <= kOrder; i++) - W_fieldnames_c[i] = W_fieldnames[i].c_str(); - plhs[2] = mxCreateStructMatrix(1, 1, kOrder+1, W_fieldnames_c); - - // Fill that structure - for (int i = 0; i <= kOrder; i++) - { - const FFSTensor &t = cond_fdr.get(Symmetry{i}); - mxArray *tmp = mxCreateDoubleMatrix(t.nrows(), t.ncols(), mxREAL); - const ConstVector &vec = t.getData(); - assert(vec.skip() == 1); - std::copy_n(vec.base(), vec.length(), mxGetPr(tmp)); - mxSetField(plhs[2], 0, ("W_" + std::to_string(i)).c_str(), tmp); - } - if (nlhs > 3) - { - - const FGSContainer &U = appwel.get_unc_ders(); - - size_t nfields = (kOrder == 1 ? 2 : (kOrder == 2 ? 6 : 12)); - const char *c_fieldnames[] = { "Uy", "Uu", "Uyy", "Uyu", "Uuu", "Uss", "Uyyy", "Uyyu", "Uyuu", "Uuuu", "Uyss", "Uuss" }; - plhs[3] = mxCreateStructMatrix(1, 1, nfields, c_fieldnames); - - copy_derivatives(plhs[3], Symmetry{1, 0, 0, 0}, U, "Uy"); - copy_derivatives(plhs[3], Symmetry{0, 1, 0, 0}, U, "Uu"); - if (kOrder >= 2) - { - copy_derivatives(plhs[3], Symmetry{2, 0, 0, 0}, U, "Uyy"); - copy_derivatives(plhs[3], Symmetry{0, 2, 0, 0}, U, "Uuu"); - copy_derivatives(plhs[3], Symmetry{1, 1, 0, 0}, U, "Uyu"); - copy_derivatives(plhs[3], Symmetry{0, 0, 0, 2}, U, "Uss"); - } - if (kOrder >= 3) - { - copy_derivatives(plhs[3], Symmetry{3, 0, 0, 0}, U, "Uyyy"); - copy_derivatives(plhs[3], Symmetry{0, 3, 0, 0}, U, "Uuuu"); - copy_derivatives(plhs[3], Symmetry{2, 1, 0, 0}, U, "Uyyu"); - copy_derivatives(plhs[3], Symmetry{1, 2, 0, 0}, U, "Uyuu"); - copy_derivatives(plhs[3], Symmetry{1, 0, 0, 2}, U, "Uyss"); - copy_derivatives(plhs[3], Symmetry{0, 1, 0, 2}, U, "Uuss"); - } - - if (nlhs > 4) - { - const FGSContainer &W = appwel.get_cond_ders(); - - size_t nfields = (kOrder == 1 ? 2 : (kOrder == 2 ? 6 : 12)); - const char *c_fieldnames[] = { "Wy", "Wu", "Wyy", "Wyu", "Wuu", "Wss", "Wyyy", "Wyyu", "Wyuu", "Wuuu", "Wyss", "Wuss" }; - plhs[4] = mxCreateStructMatrix(1, 1, nfields, c_fieldnames); - - copy_derivatives(plhs[4], Symmetry{1, 0, 0, 0}, W, "Wy"); - copy_derivatives(plhs[4], Symmetry{0, 1, 0, 0}, W, "Wu"); - if (kOrder >= 2) - { - copy_derivatives(plhs[4], Symmetry{2, 0, 0, 0}, W, "Wyy"); - copy_derivatives(plhs[4], Symmetry{0, 2, 0, 0}, W, "Wuu"); - copy_derivatives(plhs[4], Symmetry{1, 1, 0, 0}, W, "Wyu"); - copy_derivatives(plhs[4], Symmetry{0, 0, 0, 2}, W, "Wss"); - } - if (kOrder >= 3) - { - copy_derivatives(plhs[4], Symmetry{3, 0, 0, 0}, W, "Wyyy"); - copy_derivatives(plhs[4], Symmetry{0, 3, 0, 0}, W, "Wuuu"); - copy_derivatives(plhs[4], Symmetry{2, 1, 0, 0}, W, "Wyyu"); - copy_derivatives(plhs[4], Symmetry{1, 2, 0, 0}, W, "Wyuu"); - copy_derivatives(plhs[4], Symmetry{1, 0, 0, 2}, W, "Wyss"); - copy_derivatives(plhs[4], Symmetry{0, 1, 0, 2}, W, "Wuss"); - } - } - } - } - } + // Fill that structure + for (int i = 0; i <= kOrder; i++) + { + const FFSTensor &t = cond_fdr.get(Symmetry{i}); + mxArray *tmp = mxCreateDoubleMatrix(t.nrows(), t.ncols(), mxREAL); + const ConstVector &vec = t.getData(); + assert(vec.skip() == 1); + std::copy_n(vec.base(), vec.length(), mxGetPr(tmp)); + mxSetField(plhs[0], 0, ("W_" + std::to_string(i)).c_str(), tmp); + } } // end of mexFunction() } // end of extern C diff --git a/tests/optimal_policy/neo_growth_ramsey_k_order.mod b/tests/optimal_policy/neo_growth_ramsey_k_order.mod index 8e62a4399..3b994a8fd 100644 --- a/tests/optimal_policy/neo_growth_ramsey_k_order.mod +++ b/tests/optimal_policy/neo_growth_ramsey_k_order.mod @@ -31,13 +31,13 @@ end; stoch_simul(order=6, irf=0); -evaluate_planner_objective; +%evaluate_planner_objective; -[condWelfare, U_dynpp, W_dynpp, U_dyn, W_dyn] = k_order_welfare(oo_.dr, M_, options_); +[W_dynpp] = k_order_welfare(oo_.dr, M_, options_); -if condWelfare~=oo_.planner_objective_value.conditional.steady_initial_multiplier - error('Inaccurate conditional welfare with Lagrange multiplier set to its steady-state value'); -end +%if condWelfare~=oo_.planner_objective_value.conditional.steady_initial_multiplier +% error('Inaccurate conditional welfare with Lagrange multiplier set to its steady-state value'); +%end if ~exist(['neo_growth_k_order' filesep 'Output' filesep 'neo_growth_k_order_results.mat'],'file'); error('neo_growth_k_order must be run first'); end; @@ -46,20 +46,16 @@ oo = load(['neo_growth_k_order' filesep 'Output' filesep 'neo_growth_k_order_res M = load(['neo_growth_k_order' filesep 'Output' filesep 'neo_growth_k_order_results'],'M_'); options = load(['neo_growth_k_order' filesep 'Output' filesep 'neo_growth_k_order_results'],'options_'); -ind_U = strmatch('U', M.M_.endo_names,'exact'); ind_W = strmatch('W', M.M_.endo_names,'exact'); err = -1e6; for i = 1:options_.order - field_U = strcat('U_', num2str(i)); field_W = strcat('W_', num2str(i)); g_i = eval(strcat('oo.oo_.dr.g_', num2str(i))); - tmp_err = max(U_dynpp.(field_U) - g_i(ind_U, :)); - err = max(err, tmp_err); tmp_err = max(W_dynpp.(field_W) - g_i(ind_W, :)); err = max(err, tmp_err); end if err > 1e-10; - error('Inaccurate assessment of the derivatives of the felicity and the welfare functions'); + error('Inaccurate assessment of the derivatives of the welfare function'); end;