From cf829fb28e30eb1c72256afbe81209501e90c063 Mon Sep 17 00:00:00 2001 From: Normann Rion Date: Fri, 3 Dec 2021 11:07:16 +0100 Subject: [PATCH] A few fixes for k-order welfare assesment in `evaluate_planner_objective` As suggested in !1962 --- matlab/evaluate_planner_objective.m | 20 ++++++----- .../k_order_welfare/k_order_welfare.cc | 34 +++++++++---------- .../neo_growth_ramsey_k_order.mod | 2 ++ 3 files changed, 30 insertions(+), 26 deletions(-) diff --git a/matlab/evaluate_planner_objective.m b/matlab/evaluate_planner_objective.m index 65c1059b1..770738063 100644 --- a/matlab/evaluate_planner_objective.m +++ b/matlab/evaluate_planner_objective.m @@ -229,21 +229,23 @@ if options_.ramsey_policy [W] = k_order_welfare(dr,M_,options_); % Appends the welfare decision rule to the endogenous variables decision % rule - g = dr; for i=0:options_.order - eval("g.g_"+num2str(i)+" = [dr.g_"+num2str(i)+"; W.W_"+num2str(i)+"];"); + dr.(['g_' num2str(i)]) = [dr.(['g_' num2str(i)]); W.(['W_' num2str(i)])]; end % Amends the steady-state vector accordingly [U] = feval([M_.fname '.objective.static'],ys,zeros(1,exo_nbr), M_.params); ysteady = [ys(oo_.dr.order_var); U/(1-beta)]; % Generates the sequence of shocks to compute unconditional welfare - nper = 10000; - nburn = 1000; - chol_S = chol(M_.Sigma_e); - exo_simul = chol_S*randn(nper,M_.exo_nbr)'; + i_exo_var = setdiff([1:M_.exo_nbr],find(diag(M_.Sigma_e) == 0)); + nxs = length(i_exo_var); + chol_S = chol(M_.Sigma_e(i_exo_var,i_exo_var)); + exo_simul = zeros(M_.exo_nbr,options_.ramsey.nperiods); + if ~isempty(M_.Sigma_e) + exo_simul(i_exo_var,:) = chol_S*randn(nxs,options_.ramsey.nperiods); + end yhat_start = zeros(M_.endo_nbr+1,1); - [moment] = k_order_mean(options_.order, M_.nstatic, M_.npred, M_.nboth, M_.nfwrd+1, M_.exo_nbr, 1, nburn, yhat_start, exo_simul, ysteady, g); + [moment] = k_order_mean(options_.order, M_.nstatic, M_.npred, M_.nboth, M_.nfwrd+1, M_.exo_nbr, 1, options_.ramsey.drop, yhat_start, exo_simul, ysteady, dr); % Stores the result for unconditional welfare planner_objective_value.unconditional = moment(end); @@ -255,12 +257,12 @@ if options_.ramsey_policy % Conditional welfare (i) with Lagrange multipliers set to their % steady-state values yhat_start(M_.nstatic+1:M_.nstatic+M_.npred+M_.nboth) = yhat_L_SS; - [moment,sim] = k_order_mean(options_.order, M_.nstatic, M_.npred, M_.nboth, M_.nfwrd+1, M_.exo_nbr, 1, 0, yhat_start, u, ysteady, g); + [~,sim] = k_order_mean(options_.order, M_.nstatic, M_.npred, M_.nboth, M_.nfwrd+1, M_.exo_nbr, 1, 0, yhat_start, u, ysteady, dr); planner_objective_value.conditional.steady_initial_multiplier = sim(end,1); % Conditional welfare (ii) with Lagrange multipliers set to 0 yhat_start(M_.nstatic+1:M_.nstatic+M_.npred+M_.nboth) = yhat_L_0; - [moment,sim] = k_order_mean(options_.order, M_.nstatic, M_.npred, M_.nboth, M_.nfwrd+1, M_.exo_nbr, 1, nburn, yhat_start, u, ysteady, g); + [~,sim] = k_order_mean(options_.order, M_.nstatic, M_.npred, M_.nboth, M_.nfwrd+1, M_.exo_nbr, 1, 0, yhat_start, u, ysteady, dr); planner_objective_value.conditional.zero_initial_multiplier = sim(end,1); end end diff --git a/mex/sources/k_order_welfare/k_order_welfare.cc b/mex/sources/k_order_welfare/k_order_welfare.cc index 53501411e..71c57d6d4 100644 --- a/mex/sources/k_order_welfare/k_order_welfare.cc +++ b/mex/sources/k_order_welfare/k_order_welfare.cc @@ -277,25 +277,25 @@ extern "C" { ApproximationWelfare appwel(welfare, discount_factor, app.get_rule_ders(), app.get_rule_ders_s(), journal); appwel.approxAtSteady(); - const FoldDecisionRule &cond_fdr = appwel.getFoldCondWel(); - // Add possibly missing field names - for (int i = static_cast(W_fieldnames.size()); i <= kOrder; i++) + 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++) + // 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); + plhs[0] = 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[0], 0, ("W_" + std::to_string(i)).c_str(), tmp); - } + // 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 e11ff0f51..a2711b7ce 100644 --- a/tests/optimal_policy/neo_growth_ramsey_k_order.mod +++ b/tests/optimal_policy/neo_growth_ramsey_k_order.mod @@ -31,6 +31,8 @@ end; stoch_simul(order=6, irf=0); +options_.ramsey.nperiods = 10000; +options_.ramsey.drop = 1000; evaluate_planner_objective; [W_dynpp] = k_order_welfare(oo_.dr, M_, options_);