A few fixes for k-order welfare assesment in `evaluate_planner_objective`

As suggested in !1962
pac-components
Normann Rion 2021-12-03 11:07:16 +01:00 committed by Sébastien Villemot
parent 40b2565140
commit cf829fb28e
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
3 changed files with 30 additions and 26 deletions

View File

@ -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

View File

@ -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<int>(W_fieldnames.size()); i <= kOrder; i++)
const FoldDecisionRule &cond_fdr = appwel.getFoldCondWel();
// Add possibly missing field names
for (int i = static_cast<int>(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

View File

@ -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_);