Amend the k_order_welfare routine for it to return the adequate output variables
NormannR
parent
45aad05670
commit
f889a25e86
|
|
@ -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<FGSContainer>(rule_ders_arg);
|
||||
rule_ders_s = std::make_unique<FGSContainer>(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<Storage::fold>(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<FoldDecisionRule>(*unc_ders, welfare.getModel().nys(), welfare.getModel().nexog(), welfare.getModel().getSteady());
|
||||
// construct the resulting decision rule
|
||||
cond_fdr = std::make_unique<FoldDecisionRule>(*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<FGSContainer>(U);
|
||||
cond_ders = std::make_unique<FGSContainer>(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());
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -31,14 +31,9 @@ class ApproximationWelfare
|
|||
double discount_factor;
|
||||
std::unique_ptr<FGSContainer> rule_ders;
|
||||
std::unique_ptr<FGSContainer> rule_ders_s;
|
||||
std::unique_ptr<FGSContainer> unc_ders;
|
||||
std::unique_ptr<FGSContainer> cond_ders;
|
||||
std::unique_ptr<FoldDecisionRule> unc_fdr;
|
||||
std::unique_ptr<FoldDecisionRule> 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
|
||||
|
|
|
|||
|
|
@ -59,7 +59,6 @@ DynareMxArrayToString(const mxArray *mxFldp)
|
|||
)
|
||||
*/
|
||||
|
||||
std::vector<std::string> U_fieldnames;
|
||||
std::vector<std::string> 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<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++)
|
||||
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<int>(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<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++)
|
||||
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
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
|
|
|
|||
Loading…
Reference in New Issue