diff --git a/matlab/evaluate_planner_objective.m b/matlab/evaluate_planner_objective.m
index c37764ab6c8033c0e4a25297f62893204125840f..65c1059b1c5296d91b2b62365a17b12526537847 100644
--- a/matlab/evaluate_planner_objective.m
+++ b/matlab/evaluate_planner_objective.m
@@ -198,7 +198,7 @@ if options_.ramsey_policy
Wbar = U/(1-beta);
Wy = Uy*gy/(eye(nspred)-beta*Gy);
Wu = Uy*gu + beta*Wy*Gu;
-
+
if isempty(options_.qz_criterium)
options_.qz_criterium = 1+1e-6;
end
@@ -225,16 +225,43 @@ if options_.ramsey_policy
planner_objective_value.conditional.steady_initial_multiplier = W_L_SS;
planner_objective_value.conditional.zero_initial_multiplier = W_L_0;
else
- %Order k code will go here!
- if ~isempty(M_.endo_histval)
- fprintf('\nevaluate_planner_objective: order>2 conditional and unconditional welfare calculations not yet supported when an histval block is provided\n')
- else
- fprintf('\nevaluate_planner_objective: order>2 conditional welfare with initial Lagrange multipliers set to zero and unconditional welfare calculations not yet supported\n')
- planner_objective_value.conditional.steady_initial_multiplier = k_order_welfare(dr, M_, options_);
- planner_objective_value.conditional.zero_initial_multiplier = NaN;
- planner_objective_value.unconditional = NaN;
+ % Computes the welfare decision rule
+ [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)+"];");
end
- return
+ % 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)';
+ 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);
+
+ % Stores the result for unconditional welfare
+ planner_objective_value.unconditional = moment(end);
+
+ % Conditional welfare
+ % Gets initial values
+ [yhat_L_SS,yhat_L_0, u] = get_initial_state(ys,M_,dr,oo_);
+
+ % 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);
+ 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);
+ planner_objective_value.conditional.zero_initial_multiplier = sim(end,1);
end
end
elseif options_.discretionary_policy
diff --git a/mex/build/k_order_mean.am b/mex/build/k_order_mean.am
new file mode 100644
index 0000000000000000000000000000000000000000..5c38cd37c101f56d1cb927effd0247a010581662
--- /dev/null
+++ b/mex/build/k_order_mean.am
@@ -0,0 +1,14 @@
+mex_PROGRAMS = k_order_mean
+
+k_order_mean_FCFLAGS = $(AM_FCFLAGS) -I../libkordersim
+
+nodist_k_order_mean_SOURCES = \
+ mexFunction.f08
+
+k_order_mean_LDADD = ../libkordersim/libkordersim.a
+
+BUILT_SOURCES = $(nodist_k_order_mean_SOURCES)
+CLEANFILES = $(nodist_k_order_mean_SOURCES)
+
+%.f08: $(top_srcdir)/../../sources/k_order_mean/%.f08
+ $(LN_S) -f $< $@
diff --git a/mex/build/matlab/Makefile.am b/mex/build/matlab/Makefile.am
index 9609b172393ab4acdb44e9b313c5408b0051aaca..7b565fb3a07b851ca55305ebc200cfc89da512a0 100644
--- a/mex/build/matlab/Makefile.am
+++ b/mex/build/matlab/Makefile.am
@@ -4,7 +4,7 @@ SUBDIRS = mjdgges kronecker bytecode block_kalman_filter sobol perfect_foresight
# libdynare++ must come before gensylv, k_order_perturbation, dynare_simul_
if ENABLE_MEX_DYNAREPLUSPLUS
-SUBDIRS += libdynare++ gensylv libkorder dynare_simul_ k_order_perturbation k_order_welfare local_state_space_iterations libkordersim folded_to_unfolded_dr local_state_space_iteration_fortran k_order_simul
+SUBDIRS += libdynare++ gensylv libkorder dynare_simul_ k_order_perturbation k_order_welfare local_state_space_iterations libkordersim folded_to_unfolded_dr local_state_space_iteration_fortran k_order_simul k_order_mean
endif
if ENABLE_MEX_MS_SBVAR
diff --git a/mex/build/matlab/configure.ac b/mex/build/matlab/configure.ac
index 06182ea8bab8e3aade8f87472d9b6aea45f4776d..f377afc141dd5141ca1745b46deda65308f41b23 100644
--- a/mex/build/matlab/configure.ac
+++ b/mex/build/matlab/configure.ac
@@ -165,6 +165,7 @@ AC_CONFIG_FILES([Makefile
folded_to_unfolded_dr/Makefile
local_state_space_iteration_fortran/Makefile
k_order_simul/Makefile
+ k_order_mean/Makefile
perfect_foresight_problem/Makefile
num_procs/Makefile
block_trust_region/Makefile
diff --git a/mex/build/matlab/k_order_mean/Makefile.am b/mex/build/matlab/k_order_mean/Makefile.am
new file mode 100644
index 0000000000000000000000000000000000000000..58684092cff79725863fb7f3c8381ee0091992dc
--- /dev/null
+++ b/mex/build/matlab/k_order_mean/Makefile.am
@@ -0,0 +1,2 @@
+include ../mex.am
+include ../../k_order_mean.am
\ No newline at end of file
diff --git a/mex/sources/Makefile.am b/mex/sources/Makefile.am
index 214c1382fa1b94b6a28a2bec4b17d2128ad18df3..d4d924c07252022ae9905e4f8b392bcf175fe226 100644
--- a/mex/sources/Makefile.am
+++ b/mex/sources/Makefile.am
@@ -20,6 +20,7 @@ EXTRA_DIST = \
folded_to_unfolded_dr \
local_state_space_iteration_fortran \
k_order_simul \
+ k_order_mean \
gensylv \
dynare_simul_ \
perfect_foresight_problem \
diff --git a/mex/sources/k_order_mean/mexFunction.f08 b/mex/sources/k_order_mean/mexFunction.f08
new file mode 100644
index 0000000000000000000000000000000000000000..d2ecee7a5fbed63a8c1a87d3cc8ec61f06ab64d1
--- /dev/null
+++ b/mex/sources/k_order_mean/mexFunction.f08
@@ -0,0 +1,200 @@
+! Copyright © 2021 Dynare Team
+!
+! This file is part of Dynare.
+!
+! Dynare is free software: you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation, either version 3 of the License, or
+! (at your option) any later version.
+!
+! Dynare is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License
+! along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+!
+! input:
+! order the order of approximation, needs order+1 derivatives
+! nstat
+! npred
+! nboth
+! nforw
+! nexog
+! order_moment order of the moment we need to compute
+! nburn number of burn-in periods
+! yhat_start starting value for the full vector of endogenous variables minus its steady-state value
+! shocks matrix of shocks (number of exogenous variables x number of periods)
+! ysteady steady-state value for the full vector of endogenous variables
+! dr structure containing the matrices for the decision rule (g_0, g_1,…)
+! output:
+! mean estimated `order_moment`-order moment for the full vector of e
+
+subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
+ use iso_fortran_env
+ use iso_c_binding
+ use struct
+ use matlab_mex
+ use partitions
+ use simulation
+ implicit none
+
+ type(c_ptr), dimension(*), intent(in), target :: prhs
+ type(c_ptr), dimension(*), intent(out) :: plhs
+ integer(c_int), intent(in), value :: nlhs, nrhs
+ type(c_ptr) :: order_mx, nstatic_mx, npred_mx, nboth_mx, nfwrd_mx, nexog_mx, order_moment_mx, nburn_mx, &
+ yhat_start_mx, shocks_mx, ysteady_mx, dr_mx, tmp
+ type(pol), dimension(:), allocatable, target :: fdr, udr
+ integer :: order, nstatic, npred, nboth, nfwrd, exo_nbr, endo_nbr, nys, nvar, nper, nburn, order_moment
+ real(real64), dimension(:,:), allocatable :: shocks, sim
+ real(real64), dimension(:), allocatable :: dyu, mean
+ real(real64), dimension(:), pointer, contiguous :: ysteady, yhat_start
+ type(pascal_triangle) :: p
+ type(horner), dimension(:), allocatable :: h
+ integer :: i, t, d, m, n
+ character(kind=c_char, len=10) :: fieldname
+
+ order_mx = prhs(1)
+ nstatic_mx = prhs(2)
+ npred_mx = prhs(3)
+ nboth_mx = prhs(4)
+ nfwrd_mx = prhs(5)
+ nexog_mx = prhs(6)
+ order_moment_mx = prhs(7)
+ nburn_mx = prhs(8)
+ yhat_start_mx = prhs(9)
+ shocks_mx = prhs(10)
+ ysteady_mx = prhs(11)
+ dr_mx = prhs(12)
+
+ ! Checking the consistence and validity of input arguments
+ if (nrhs /= 12) then
+ call mexErrMsgTxt("Must have exactly 12 inputs")
+ end if
+ if (.not. (mxIsScalar(order_mx)) .and. mxIsNumeric(order_mx)) then
+ call mexErrMsgTxt("1st argument (order) should be a numeric scalar")
+ end if
+ if (.not. (mxIsScalar(nstatic_mx)) .and. mxIsNumeric(nstatic_mx)) then
+ call mexErrMsgTxt("2nd argument (nstat) should be a numeric scalar")
+ end if
+ if (.not. (mxIsScalar(npred_mx)) .and. mxIsNumeric(npred_mx)) then
+ call mexErrMsgTxt("3rd argument (npred) should be a numeric scalar")
+ end if
+ if (.not. (mxIsScalar(nboth_mx)) .and. mxIsNumeric(nboth_mx)) then
+ call mexErrMsgTxt("4th argument (nboth) should be a numeric scalar")
+ end if
+ if (.not. (mxIsScalar(nfwrd_mx)) .and. mxIsNumeric(nfwrd_mx)) then
+ call mexErrMsgTxt("5th argument (nforw) should be a numeric scalar")
+ end if
+ if (.not. (mxIsScalar(nexog_mx)) .and. mxIsNumeric(nexog_mx)) then
+ call mexErrMsgTxt("6th argument (nexog) should be a numeric scalar")
+ end if
+ if (.not. (mxIsScalar(order_moment_mx)) .and. mxIsNumeric(order_moment_mx)) then
+ call mexErrMsgTxt("7th argument (order_moment) should be a numeric scalar")
+ end if
+ if (.not. (mxIsScalar(nburn_mx)) .and. mxIsNumeric(nburn_mx)) then
+ call mexErrMsgTxt("8th argument (nburn) should be a numeric scalar")
+ end if
+ if (.not. (mxIsDouble(yhat_start_mx) .and. (mxGetM(yhat_start_mx) == 1 .or. mxGetN(yhat_start_mx) == 1))) then
+ call mexErrMsgTxt("9th argument (yhat_start) should be a real vector")
+ end if
+ if (.not. (mxIsDouble(shocks_mx))) then
+ call mexErrMsgTxt("10th argument (shocks) should be a real matrix")
+ end if
+ if (.not. (mxIsDouble(ysteady_mx) .and. (mxGetM(ysteady_mx) == 1 .or. mxGetN(ysteady_mx) == 1))) then
+ call mexErrMsgTxt("11th argument (ysteady) should be a real vector")
+ end if
+ if (.not. mxIsStruct(dr_mx)) then
+ call mexErrMsgTxt("12th argument (dr) should be a struct")
+ end if
+
+ ! Converting inputs in Fortran format
+ order = int(mxGetScalar(order_mx))
+ nstatic = int(mxGetScalar(nstatic_mx))
+ npred = int(mxGetScalar(npred_mx))
+ nboth = int(mxGetScalar(nboth_mx))
+ nfwrd = int(mxGetScalar(nfwrd_mx))
+ exo_nbr = int(mxGetScalar(nexog_mx))
+ endo_nbr = nstatic+npred+nboth+nfwrd
+ nys = npred+nboth
+ nvar = nys+exo_nbr
+ nburn = int(mxGetScalar(nburn_mx))
+ order_moment = int(mxGetScalar(order_moment_mx))
+
+ if (endo_nbr /= int(mxGetM(yhat_start_mx))) then
+ call mexErrMsgTxt("yhat_start should have nstat+npred+nboth+nforw rows")
+ end if
+ yhat_start => mxGetPr(yhat_start_mx)
+
+ if (exo_nbr /= int(mxGetM(shocks_mx))) then
+ call mexErrMsgTxt("shocks should have nexog rows")
+ end if
+ nper = int(mxGetN(shocks_mx))
+ allocate(shocks(exo_nbr,nper))
+ shocks = reshape(mxGetPr(shocks_mx),[exo_nbr,nper])
+
+ if (.not. (int(mxGetM(ysteady_mx)) == endo_nbr)) then
+ call mexErrMsgTxt("ysteady should have nstat+npred+nboth+nforw rows")
+ end if
+ ysteady => mxGetPr(ysteady_mx)
+
+ allocate(h(0:order), fdr(0:order), udr(0:order))
+ do i = 0, order
+ write (fieldname, '(a2, i1)') "g_", i
+ tmp = mxGetField(dr_mx, 1_mwIndex, trim(fieldname))
+ if (.not. (c_associated(tmp) .and. mxIsDouble(tmp))) then
+ call mexErrMsgTxt(trim(fieldname)//" is not allocated in dr")
+ end if
+ m = int(mxGetM(tmp))
+ n = int(mxGetN(tmp))
+ allocate(fdr(i)%g(m,n), udr(i)%g(endo_nbr, nvar**i), h(i)%c(endo_nbr, nvar**i))
+ fdr(i)%g(1:m,1:n) = reshape(mxGetPr(tmp), [m,n])
+ end do
+
+ udr(0)%g = fdr(0)%g
+ udr(1)%g = fdr(1)%g
+ if (order > 1) then
+ ! Compute the useful binomial coefficients from Pascal's triangle
+ p = pascal_triangle(nvar+order-1)
+ block
+ type(uf_matching), dimension(2:order) :: matching
+ ! Pinpointing the corresponding offsets between folded and unfolded tensors
+ do d=2,order
+ allocate(matching(d)%folded(nvar**d))
+ call fill_folded_indices(matching(d)%folded, nvar, d, p)
+ udr(d)%g = fdr(d)%g(:,matching(d)%folded)
+ end do
+ end block
+ end if
+
+ allocate(dyu(nvar), mean(endo_nbr), sim(endo_nbr,nper))
+ ! Getting the predetermined part of the endogenous variable vector
+ dyu(1:nys) = yhat_start(nstatic+1:nstatic+nys)
+ dyu(nys+1:) = shocks(:,1)
+ ! Using the Horner algorithm to evaluate the decision rule at the chosen dyu
+ call eval(h, dyu, udr, endo_nbr, nvar, order)
+ sim(:,1) = h(0)%c(:,1) + ysteady
+ mean = 0.
+
+ ! Carrying out the simulation
+ do t=2,nper
+ dyu(1:nys) = h(0)%c(nstatic+1:nstatic+nys,1)
+ dyu(nys+1:) = shocks(:,t)
+ call eval(h, dyu, udr, endo_nbr, nvar, order)
+ sim(:,t) = h(0)%c(:,1) + ysteady
+ if (t > nburn) then
+ mean = mean + sim(:,t)**order_moment
+ end if
+ end do
+ ! scaling the mean with the number of non-burn-in periods
+ mean = mean/(nper-nburn)
+
+ ! Generating output
+ plhs(1) = mxCreateDoubleMatrix(int(endo_nbr, mwSize), 1_mwSize, mxREAL)
+ mxGetPr(plhs(1)) = mean
+ plhs(2) = mxCreateDoubleMatrix(int(endo_nbr, mwSize), int(nper, mwSize), mxREAL)
+ mxGetPr(plhs(2)) = reshape(sim, [size(sim)])
+
+
+end subroutine mexFunction
diff --git a/mex/sources/k_order_welfare/approximation_welfare.cc b/mex/sources/k_order_welfare/approximation_welfare.cc
index 5e3b527f9c6fd6c8e65df24e847448e026b6490c..ccd82854cf1bfc39b34fcf7aaf6782eaec71d8b3 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<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());
- }
-}
+}
\ 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 c2d2d166f9e8c8d67933fac868772a367ed5a62b..9c22eae259edc15b4350fb4afcd5ad54543111f2 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<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
diff --git a/mex/sources/k_order_welfare/k_order_welfare.cc b/mex/sources/k_order_welfare/k_order_welfare.cc
index 266d133383553a75bbe70e1bd095e3de4acc07f7..53501411eccfe0fea199cc4207d72462d77a72cd 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<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]));
-
- 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");
- }
- }
- }
- }
- }
+ 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);
+
+ // 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/mex/sources/local_state_space_iteration_fortran/mexFunction.f08 b/mex/sources/local_state_space_iteration_fortran/mexFunction.f08
index ae3662aa2094192d3d979eff579ca1f6fa6929f6..9d4d32af33f104ef7b5d3a7760bd00383889c598 100644
--- a/mex/sources/local_state_space_iteration_fortran/mexFunction.f08
+++ b/mex/sources/local_state_space_iteration_fortran/mexFunction.f08
@@ -17,7 +17,7 @@
! input:
! yhat values of endogenous variables
-! epsilon values of the exgogenous shock
+! epsilon values of the exogenous shock
! dr struct containing the folded tensors g_0, g_1, ...
! M struct containing the model features
! options struct containing the model options
diff --git a/tests/optimal_policy/neo_growth_ramsey_k_order.mod b/tests/optimal_policy/neo_growth_ramsey_k_order.mod
index 8e62a43992043180633e059aa9c30c4d2a6c2bf4..e11ff0f519361f200592095a1fca3f5e4eb6abc8 100644
--- a/tests/optimal_policy/neo_growth_ramsey_k_order.mod
+++ b/tests/optimal_policy/neo_growth_ramsey_k_order.mod
@@ -33,11 +33,8 @@ stoch_simul(order=6, irf=0);
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 ~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 +43,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;