diff --git a/dynare++/kord/approximation.cc b/dynare++/kord/approximation.cc
index 0d082edb367406e17976669a88dc656fe0d0329c..c10be1696c82063f8100d4fa9380badf75924495 100644
--- a/dynare++/kord/approximation.cc
+++ b/dynare++/kord/approximation.cc
@@ -1,6 +1,6 @@
/*
* Copyright © 2005 Ondra Kamenik
- * Copyright © 2019 Dynare Team
+ * Copyright © 2019-2021 Dynare Team
*
* This file is part of Dynare.
*
@@ -240,11 +240,14 @@ void
Approximation::saveRuleDerivs(const FGSContainer &g)
{
rule_ders = std::make_unique<FGSContainer>(g);
+ rule_ders_s = std::make_unique<FGSContainer>(4);
rule_ders_ss = std::make_unique<FGSContainer>(4);
for (auto &run : *rule_ders)
{
- auto ten = std::make_unique<FGSTensor>(ypart.nstat+ypart.npred, ypart.nyss(), *(run.second));
- rule_ders_ss->insert(std::move(ten));
+ auto ten_s = std::make_unique<FGSTensor>(ypart.nstat, ypart.nys(), *(run.second));
+ rule_ders_s->insert(std::move(ten_s));
+ auto ten_ss = std::make_unique<FGSTensor>(ypart.nstat+ypart.npred, ypart.nyss(), *(run.second));
+ rule_ders_ss->insert(std::move(ten_ss));
}
}
diff --git a/dynare++/kord/approximation.hh b/dynare++/kord/approximation.hh
index fdf0bd5fa0883568ac410491e5d094bc9566248b..42b31290193e0928a43cba6ed478ddc1423b1113 100644
--- a/dynare++/kord/approximation.hh
+++ b/dynare++/kord/approximation.hh
@@ -1,6 +1,6 @@
/*
* Copyright © 2005 Ondra Kamenik
- * Copyright © 2019 Dynare Team
+ * Copyright © 2019-2021 Dynare Team
*
* This file is part of Dynare.
*
@@ -122,6 +122,7 @@ class Approximation
DynamicModel &model;
Journal &journal;
std::unique_ptr<FGSContainer> rule_ders;
+ std::unique_ptr<FGSContainer> rule_ders_s;
std::unique_ptr<FGSContainer> rule_ders_ss;
std::unique_ptr<FoldDecisionRule> fdr;
std::unique_ptr<UnfoldDecisionRule> udr;
@@ -156,6 +157,11 @@ public:
return *rule_ders;
}
const FGSContainer &
+ get_rule_ders_s() const
+ {
+ return *rule_ders_s;
+ }
+ const FGSContainer &
get_rule_ders_ss() const
{
return *rule_ders_ss;
diff --git a/dynare++/kord/decision_rule.hh b/dynare++/kord/decision_rule.hh
index 23b4636589aabb278333f26bfa50be72a57d0749..2736ab4cc4cc092d120a833faee9b71fc5a97919 100644
--- a/dynare++/kord/decision_rule.hh
+++ b/dynare++/kord/decision_rule.hh
@@ -1,6 +1,6 @@
/*
* Copyright © 2004 Ondra Kamenik
- * Copyright © 2019 Dynare Team
+ * Copyright © 2019-2021 Dynare Team
*
* This file is part of Dynare.
*
@@ -126,6 +126,7 @@ class DecisionRuleImpl : public ctraits<t>::Tpol, public DecisionRule
protected:
using _Tpol = typename ctraits<t>::Tpol;
using _Tg = typename ctraits<t>::Tg;
+ using _TW = typename ctraits<t>::TW;
using _Ttensor = typename ctraits<t>::Ttensor;
using _Ttensym = typename ctraits<t>::Ttensym;
const Vector ysteady;
@@ -149,6 +150,12 @@ public:
{
fillTensors(g, sigma);
}
+ DecisionRuleImpl(const _TW &W, int nys, int nuu,
+ const ConstVector &ys)
+ : ctraits<t>::Tpol(1, nys+nuu), ysteady(ys), nu(nuu)
+ {
+ fillTensors(W, nys);
+ }
DecisionRuleImpl(const DecisionRuleImpl<t> &dr, const ConstVector &fixpoint)
: ctraits<t>::Tpol(dr.ypart.ny(), dr.ypart.nys()+dr.nu),
ysteady(fixpoint), ypart(dr.ypart), nu(dr.nu)
@@ -179,6 +186,7 @@ public:
}
protected:
void fillTensors(const _Tg &g, double sigma);
+ void fillTensors(const _TW &W, int nys);
void centralize(const DecisionRuleImpl &dr);
public:
void eval(emethod em, Vector &out, const ConstVector &v) const override;
@@ -265,6 +273,50 @@ DecisionRuleImpl<t>::fillTensors(const _Tg &g, double sigma)
}
}
+template<Storage t>
+void
+DecisionRuleImpl<t>::fillTensors(const _TW &W, int nys)
+{
+ IntSequence tns{nys, nu};
+ int dfact = 1;
+ for (int d = 0; d <= W.getMaxDim(); d++, dfact *= d)
+ {
+ auto W_yud = std::make_unique<_Ttensym>(1, nys+nu, d);
+ W_yud->zeros();
+
+ // fill tensor of ‘g_yud’ of dimension ‘d’
+ /* Here we have to fill the tensor [g_(yu)ᵈ]. So we go through all pairs
+ (i,j) such that i+j=d, and through all k from zero up to maximal
+ dimension minus d. In this way we go through all symmetries of
+ [g_yⁱuʲσᵏ] which will be added to [g_(yu)ᵈ].
+
+ Note that at the beginning, ‘dfact’ is a factorial of ‘d’. We
+ calculate ‘kfact’ is equal to k!. As indicated in
+ DecisionRuleImpl::fillTensors(), the added tensor is thus multiplied
+ with 1/(d!k!)·σᵏ. */
+
+ for (int i = 0; i <= d; i++)
+ {
+ int j = d-i;
+ int kfact = 1;
+ _Ttensor tmp(1, TensorDimens(Symmetry{i, j}, tns));
+ tmp.zeros();
+ for (int k = 0; k+d <= W.getMaxDim(); k++, kfact *= k)
+ {
+ Symmetry sym{i, j, 0, k};
+ if (W.check(sym))
+ {
+ double mult = 1.0/dfact/kfact;
+ tmp.add(mult, W.get(sym));
+ }
+ }
+ W_yud->addSubTensor(tmp);
+ }
+
+ this->insert(std::move(W_yud));
+ }
+}
+
/* The centralization is straightforward. We suppose here that the object’s
steady state is the fix point ỹ. It is clear that the new derivatives
[g~_(yu)ⁱ] will be equal to the derivatives of the original decision rule
@@ -457,6 +509,11 @@ public:
: DecisionRuleImpl<Storage::fold>(g, yp, nuu, ys, sigma)
{
}
+ FoldDecisionRule(const ctraits<Storage::fold>::TW &W, int nys, int nuu,
+ const ConstVector &ys)
+ : DecisionRuleImpl<Storage::fold>(W, nys, nuu, ys)
+ {
+ }
FoldDecisionRule(const DecisionRuleImpl<Storage::fold> &dr, const ConstVector &fixpoint)
: DecisionRuleImpl<Storage::fold>(dr, fixpoint)
{
diff --git a/dynare++/kord/korder.hh b/dynare++/kord/korder.hh
index 54a098b0795578e4813e9fdfa397a66b88513f18..09d10bc826e8a23206032e2fc956fdded76937b4 100644
--- a/dynare++/kord/korder.hh
+++ b/dynare++/kord/korder.hh
@@ -1,6 +1,6 @@
/*
* Copyright © 2004 Ondra Kamenik
- * Copyright © 2019 Dynare Team
+ * Copyright © 2019-2021 Dynare Team
*
* This file is part of Dynare.
*
@@ -82,7 +82,11 @@ public:
using Tgs = std::conditional_t<type == Storage::fold, FGSContainer, UGSContainer>;
using Tgss = std::conditional_t<type == Storage::fold, FGSContainer, UGSContainer>;
using TG = std::conditional_t<type == Storage::fold, FGSContainer, UGSContainer>;
+ using TU = std::conditional_t<type == Storage::fold, FGSContainer, UGSContainer>;
+ using TW = std::conditional_t<type == Storage::fold, FGSContainer, UGSContainer>;
+ using TWrond = std::conditional_t<type == Storage::fold, FGSContainer, UGSContainer>;
using TZstack = std::conditional_t<type == Storage::fold, FoldedZContainer, UnfoldedZContainer>;
+ using TXstack = std::conditional_t<type == Storage::fold, FoldedXContainer, UnfoldedXContainer>;
using TGstack = std::conditional_t<type == Storage::fold, FoldedGContainer, UnfoldedGContainer>;
using Tm = std::conditional_t<type == Storage::fold, FNormalMoments, UNormalMoments>;
using Tpol = std::conditional_t<type == Storage::fold, FTensorPolynomial, UTensorPolynomial>;
@@ -113,6 +117,7 @@ struct PartitionY
const int npred;
const int nboth;
const int nforw;
+ PartitionY() : PartitionY(0,0,0,0) {}
PartitionY(int num_stat, int num_pred,
int num_both, int num_forw)
: nstat(num_stat), npred(num_pred),
@@ -322,6 +327,11 @@ public:
{
return _ug;
}
+ const FGSContainer &
+ getFoldDersS() const
+ {
+ return _fgs;
+ }
static bool
is_even(int i)
{
diff --git a/dynare++/kord/korder_stoch.hh b/dynare++/kord/korder_stoch.hh
index cbf55b1224b63d074b9c14d328f60a153634b0e8..79f496b366fbde8d628ad464646f9c17a532a3aa 100644
--- a/dynare++/kord/korder_stoch.hh
+++ b/dynare++/kord/korder_stoch.hh
@@ -1,6 +1,6 @@
/*
* Copyright © 2005 Ondra Kamenik
- * Copyright © 2019 Dynare Team
+ * Copyright © 2019-2021 Dynare Team
*
* This file is part of Dynare.
*
@@ -458,6 +458,11 @@ public:
{
return _ug;
}
+ const FGSContainer &
+ getFoldDersS() const
+ {
+ return _fgs;
+ }
protected:
template<Storage t>
std::unique_ptr<typename ctraits<t>::Ttensor> faaDiBrunoZ(const Symmetry &sym) const;
diff --git a/dynare++/tl/cc/stack_container.hh b/dynare++/tl/cc/stack_container.hh
index 68a61dbdb79857206f00314d20994177bd28202d..dd0c32752876e32f684b6c483995b740b1798c0a 100644
--- a/dynare++/tl/cc/stack_container.hh
+++ b/dynare++/tl/cc/stack_container.hh
@@ -337,6 +337,65 @@ protected:
UGSTensor &out, std::mutex &mut) const;
};
+/* Here is the specialization of the StackContainer. We implement
+ here the x needed in DSGE context for welfare assessment. We implement getType() and define a constructor feeding the data and sizes.
+
+ It depends on four variables U(y,u,u',σ), the variable u' being introduced to enable additions with 4-variable tensors*/
+
+template<class _Ttype>
+class XContainer : public StackContainer<_Ttype>
+{
+public:
+ using _Tparent = StackContainer<_Ttype>;
+ using _Stype = StackContainerInterface<_Ttype>;
+ using _Ctype = typename _Tparent::_Ctype;
+ using itype = typename _Tparent::itype;
+ XContainer(const _Ctype *g, int ng)
+ : _Tparent(1, 1)
+ {
+ _Tparent::stack_sizes = { ng };
+ _Tparent::conts[0] = g;
+ _Tparent::calculateOffsets();
+ }
+
+ /* Here we say, what happens if we derive z. recall the top of the
+ file, how z looks, and code is clear. */
+
+ itype
+ getType(int i, const Symmetry &s) const override
+ {
+ if (i==0)
+ if (s[2] > 0)
+ return itype::zero;
+ else
+ return itype::matrix;
+
+ TL_RAISE("Wrong stack index in XContainer::getType");
+ }
+
+};
+
+class FoldedXContainer : public XContainer<FGSTensor>,
+ public FoldedStackContainer
+{
+public:
+ using _Ctype = TensorContainer<FGSTensor>;
+ FoldedXContainer(const _Ctype *g, int ng)
+ : XContainer<FGSTensor>(g, ng)
+ {
+ }
+};
+
+class UnfoldedXContainer : public XContainer<UGSTensor>,
+ public UnfoldedStackContainer
+{
+public:
+ using _Ctype = TensorContainer<UGSTensor>;
+ UnfoldedXContainer(const _Ctype *g, int ng)
+ : XContainer<UGSTensor>(g, ng)
+ {
+ }
+};
/* Here is the specialization of the StackContainer. We implement
here the z needed in DSGE context. We implement getType() and define
a constructor feeding the data and sizes.
diff --git a/matlab/evaluate_planner_objective.m b/matlab/evaluate_planner_objective.m
index 57d7df89cd464f09955c4fe3dc2435dd651db044..d912c9cd171b7cf389b6f646c59515fc28be86b4 100644
--- a/matlab/evaluate_planner_objective.m
+++ b/matlab/evaluate_planner_objective.m
@@ -169,8 +169,8 @@ if options_.ramsey_policy
planner_objective_value(2) = W;
else
%Order k code will go here!
- fprintf('\nevaluate_planner_objective: order>2 not yet supported\n')
- planner_objective_value(1) = NaN;
+ fprintf('\nevaluate_planner_objective: order>2 unconditional welfare calculation not yet supported\n')
+ planner_objective_value(1) = k_order_welfare(dr, M_, options_);
planner_objective_value(2) = NaN;
return
end
diff --git a/matlab/mex/k_order_welfare.m b/matlab/mex/k_order_welfare.m
new file mode 100644
index 0000000000000000000000000000000000000000..95a1c3447c9e0373594c9744a2787c5ec013ab87
--- /dev/null
+++ b/matlab/mex/k_order_welfare.m
@@ -0,0 +1,60 @@
+% [unc_welfare, U_dynpp_derivs, W_dynpp_derivs, U_dyn_derivs, W_dyn_derivs] = k_order_welfare(dr, DynareModel, DynareOptions)
+% computes a k-th order approximation of welfare
+%
+% INPUTS
+% dr: struct describing the reduced form solution of the model.
+% DynareModel: struct jobs's parameters
+% DynareOptions: struct job's options
+%
+% OUTPUTS
+%
+% unc_welfare double approximation of conditional welfare from the non-stochastic steady state and allowing stochastic shocks
+%
+% U_dynpp_derivs struct Derivatives of the felicity function in Dynare++ format.
+% The tensors are folded and the Taylor coefficients (1/n!) are
+% included.
+% Fieldnames are U_0, U_1, U_2, …
+% W_dynpp_derivs struct Derivatives of the welfare function in Dynare++ format.
+% The tensors are folded and the Taylor coefficients (1/n!) are
+% included.
+% Fieldnames are W_0, W_1, W_2, …
+% U_dyn_derivs struct Derivatives of the felicity function in Dynare format, up to third order.
+% Matrices are dense and unfolded. The Taylor coefficients (1/2
+% and 1/6) aren't included.
+% The derivatives w.r.t. different categories of variables
+% are separated.
+% Fieldnames are:
+% + Uy, Uu
+% + if order ≥ 2: Uyy, Uyu, Uuu, Uss
+% + if order ≥ 3: Uyyy, Uyyu, Uyuu, Uuuu, Uyss, Uuss
+%
+% W_dyn_derivs struct Derivatives of the welfare function in Dynare format, up to third order.
+% Matrices are dense and unfolded. The Taylor coefficients (1/2
+% and 1/6) aren't included.
+% The derivatives w.r.t. different categories of variables
+% are separated.
+% Fieldnames are:
+% + Wy, Wu
+% + if order ≥ 2: Wyy, Wyu, Wuu, Wss
+% + if order ≥ 3: Wyyy, Wyyu, Wyuu, Wuuu, Wyss, Wuss
+%
+% k_order_welfare is a compiled MEX function. Its source code is in
+% dynare/mex/sources/k_order_welfare/k_order_welfare.cc and it uses code provided by
+% dynare++
+
+% Copyright (C) 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/>.
diff --git a/mex/build/k_order_perturbation.am b/mex/build/k_order_perturbation.am
index 5492a6e9133b9dbfe82b26403c7328c43d762fdf..9aed1edf793a52c39d135c9aafb742c56575eb7d 100644
--- a/mex/build/k_order_perturbation.am
+++ b/mex/build/k_order_perturbation.am
@@ -7,13 +7,9 @@ k_order_perturbation_CPPFLAGS = $(AM_CPPFLAGS) -I$(top_srcdir)/../../../dynare++
k_order_perturbation_CXXFLAGS = $(AM_CXXFLAGS) $(THREAD_CXXFLAGS)
k_order_perturbation_LDFLAGS = $(AM_LDFLAGS) $(LDFLAGS_MATIO)
-k_order_perturbation_LDADD = ../libdynare++/libdynare++.a $(LIBADD_DLOPEN) $(LIBADD_MATIO)
+k_order_perturbation_LDADD = ../libdynare++/libdynare++.a ../libkorder/libkorder.a $(LIBADD_DLOPEN) $(LIBADD_MATIO)
-nodist_k_order_perturbation_SOURCES = \
- k_order_perturbation.cc \
- k_ord_dynare.cc \
- dynamic_dll.cc \
- dynamic_m.cc
+nodist_k_order_perturbation_SOURCES = k_order_perturbation.cc
BUILT_SOURCES = $(nodist_k_order_perturbation_SOURCES)
CLEANFILES = $(nodist_k_order_perturbation_SOURCES)
diff --git a/mex/build/k_order_welfare.am b/mex/build/k_order_welfare.am
new file mode 100644
index 0000000000000000000000000000000000000000..62da9fe2b5bc4daca39327eb70da52dc6a32ed8a
--- /dev/null
+++ b/mex/build/k_order_welfare.am
@@ -0,0 +1,22 @@
+mex_PROGRAMS = k_order_welfare
+
+TOPDIR = $(top_srcdir)/../../sources/k_order_welfare
+
+k_order_welfare_CPPFLAGS = $(AM_CPPFLAGS) -I$(TOPDIR) -I$(TOPDIR)/../k_order_perturbation -I$(top_srcdir)/../../../dynare++/tl/cc -I$(top_srcdir)/../../../dynare++/kord -I$(top_srcdir)/../../../dynare++/src -I$(top_srcdir)/../../../dynare++/sylv/cc -I$(top_srcdir)/../../../dynare++/utils/cc $(CPPFLAGS_MATIO)
+
+k_order_welfare_CXXFLAGS = $(AM_CXXFLAGS) $(THREAD_CXXFLAGS)
+
+k_order_welfare_LDFLAGS = $(AM_LDFLAGS) $(LDFLAGS_MATIO)
+k_order_welfare_LDADD = ../libdynare++/libdynare++.a ../libkorder/libkorder.a $(LIBADD_DLOPEN) $(LIBADD_MATIO)
+
+nodist_k_order_welfare_SOURCES = \
+ k_order_welfare.cc \
+ approximation_welfare.cc \
+ k_ord_objective.cc \
+ objective_m.cc
+
+BUILT_SOURCES = $(nodist_k_order_welfare_SOURCES)
+CLEANFILES = $(nodist_k_order_welfare_SOURCES)
+
+%.cc: $(TOPDIR)/%.cc
+ $(LN_S) -f $< $@
diff --git a/mex/build/libkorder.am b/mex/build/libkorder.am
new file mode 100644
index 0000000000000000000000000000000000000000..59f023c84808ccac40cba6e6a0f3ec79b9dde6e4
--- /dev/null
+++ b/mex/build/libkorder.am
@@ -0,0 +1,18 @@
+noinst_LIBRARIES = libkorder.a
+
+TOPDIR = $(top_srcdir)/../../sources/k_order_perturbation
+
+libkorder_a_CPPFLAGS = $(AM_CPPFLAGS) -I$(TOPDIR) -I$(top_srcdir)/../../../dynare++/tl/cc -I$(top_srcdir)/../../../dynare++/sylv/cc -I$(top_srcdir)/../../../dynare++/src -I$(top_srcdir)/../../../dynare++/kord -I$(top_srcdir)/../../../dynare++/utils/cc $(CPPFLAGS_MATIO)
+
+libkorder_a_CXXFLAGS = $(AM_CXXFLAGS) $(THREAD_CXXFLAGS)
+
+nodist_libkorder_a_SOURCES = \
+ k_ord_dynare.cc \
+ dynamic_dll.cc \
+ dynamic_m.cc
+
+BUILT_SOURCES = $(nodist_libkorder_a_SOURCES)
+CLEANFILES = $(nodist_libkorder_a_SOURCES)
+
+%.cc: $(TOPDIR)/%.cc
+ $(LN_S) -f $< $@
diff --git a/mex/build/matlab/Makefile.am b/mex/build/matlab/Makefile.am
index e84baee6ed388fe01510b7a2d7a63d0e045b98b4..58f376cdc5e0614ec03d96f8e39c2909a302cd38 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 k_order_perturbation dynare_simul_ local_state_space_iterations
+SUBDIRS += libdynare++ gensylv libkorder dynare_simul_ k_order_perturbation k_order_welfare local_state_space_iterations
endif
if ENABLE_MEX_MS_SBVAR
diff --git a/mex/build/matlab/configure.ac b/mex/build/matlab/configure.ac
index b6ab2178e64ce458483f95005c1e02b9b714fa06..e81dd217e00fa6dd30364660137ee163f41c9556 100644
--- a/mex/build/matlab/configure.ac
+++ b/mex/build/matlab/configure.ac
@@ -153,7 +153,9 @@ AC_CONFIG_FILES([Makefile
bytecode/Makefile
libdynare++/Makefile
gensylv/Makefile
+ libkorder/Makefile
k_order_perturbation/Makefile
+ k_order_welfare/Makefile
dynare_simul_/Makefile
kalman_steady_state/Makefile
ms_sbvar/Makefile
diff --git a/mex/build/matlab/k_order_welfare/Makefile.am b/mex/build/matlab/k_order_welfare/Makefile.am
new file mode 100644
index 0000000000000000000000000000000000000000..ea0c407891403b3d8892623e8ddd77246c5ef722
--- /dev/null
+++ b/mex/build/matlab/k_order_welfare/Makefile.am
@@ -0,0 +1,2 @@
+include ../mex.am
+include ../../k_order_welfare.am
diff --git a/mex/build/matlab/libkorder/Makefile.am b/mex/build/matlab/libkorder/Makefile.am
new file mode 100644
index 0000000000000000000000000000000000000000..bc7168b3315cace306e576625d5a37cbb9d8c192
--- /dev/null
+++ b/mex/build/matlab/libkorder/Makefile.am
@@ -0,0 +1,2 @@
+include ../mex.am
+include ../../libkorder.am
diff --git a/mex/build/octave/Makefile.am b/mex/build/octave/Makefile.am
index 6726be71e928d36dc9cac5cb8deebed59ebf1b60..a669cf8393c3bd88866dc73c79f7f0dae6084df3 100644
--- a/mex/build/octave/Makefile.am
+++ b/mex/build/octave/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 k_order_perturbation dynare_simul_ local_state_space_iterations
+SUBDIRS += libdynare++ gensylv libkorder dynare_simul_ k_order_perturbation k_order_welfare local_state_space_iterations
endif
if ENABLE_MEX_MS_SBVAR
diff --git a/mex/build/octave/configure.ac b/mex/build/octave/configure.ac
index 7df4785c35c60a99f98c11f7cd930adc7fc957c5..4ac67e6668775ad0df7dc37303e14465a3213a33 100644
--- a/mex/build/octave/configure.ac
+++ b/mex/build/octave/configure.ac
@@ -135,7 +135,9 @@ AC_CONFIG_FILES([Makefile
bytecode/Makefile
libdynare++/Makefile
gensylv/Makefile
+ libkorder/Makefile
k_order_perturbation/Makefile
+ k_order_welfare/Makefile
dynare_simul_/Makefile
kalman_steady_state/Makefile
ms_sbvar/Makefile
diff --git a/mex/build/octave/k_order_welfare/Makefile.am b/mex/build/octave/k_order_welfare/Makefile.am
new file mode 100644
index 0000000000000000000000000000000000000000..4e3f52f20c373c8c98d5c08be099a7330bcbd113
--- /dev/null
+++ b/mex/build/octave/k_order_welfare/Makefile.am
@@ -0,0 +1,3 @@
+EXEEXT = .mex
+include ../mex.am
+include ../../k_order_welfare.am
diff --git a/mex/build/octave/libkorder/Makefile.am b/mex/build/octave/libkorder/Makefile.am
new file mode 100644
index 0000000000000000000000000000000000000000..6f0aef6218a710c80af506f296c27e641c820906
--- /dev/null
+++ b/mex/build/octave/libkorder/Makefile.am
@@ -0,0 +1,3 @@
+EXEEXT = .mex
+include ../mex.am
+include ../../libkorder.am
diff --git a/mex/sources/Makefile.am b/mex/sources/Makefile.am
index 3fca510decf6c3826ff2a9ea512e84e32d251b9a..1840bbf6a39724ae6782d6bd76c4858c6aad79e4 100644
--- a/mex/sources/Makefile.am
+++ b/mex/sources/Makefile.am
@@ -10,6 +10,7 @@ EXTRA_DIST = \
kronecker \
bytecode \
k_order_perturbation \
+ k_order_welfare \
kalman_steady_state \
ms-sbvar \
block_kalman_filter \
diff --git a/mex/sources/k_order_perturbation/k_ord_dynare.hh b/mex/sources/k_order_perturbation/k_ord_dynare.hh
index 4f82d46c55b776be0eb9eb023281e09a69418bf2..9cccbef7ff3afcb6440e439aa2473602ee59abf8 100644
--- a/mex/sources/k_order_perturbation/k_ord_dynare.hh
+++ b/mex/sources/k_order_perturbation/k_ord_dynare.hh
@@ -138,10 +138,30 @@ public:
return nExog;
}
int
+ ny() const
+ {
+ return nStat+nBoth+nPred+nForw;
+ }
+ int
+ nys() const
+ {
+ return nBoth+nPred;
+ }
+ int
order() const override
{
return nOrder;
}
+ const std::vector<int> &
+ getDynppToDyn() const
+ {
+ return dynppToDyn;
+ }
+ const std::vector<int> &
+ getDynToDynpp() const
+ {
+ return dynToDynpp;
+ }
const NameList &
getAllEndoNames() const override
{
diff --git a/mex/sources/k_order_welfare/approximation_welfare.cc b/mex/sources/k_order_welfare/approximation_welfare.cc
new file mode 100644
index 0000000000000000000000000000000000000000..5e3b527f9c6fd6c8e65df24e847448e026b6490c
--- /dev/null
+++ b/mex/sources/k_order_welfare/approximation_welfare.cc
@@ -0,0 +1,96 @@
+/*
+ * Copyright © 2005 Ondra Kamenik
+ * Copyright © 2019-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/>.
+ */
+
+#include <utility>
+
+#include "kord_exception.hh"
+#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}
+{
+ rule_ders = std::make_unique<FGSContainer>(rule_ders_arg);
+ rule_ders_s = std::make_unique<FGSContainer>(rule_ders_s_arg);
+}
+
+/* This methods assumes that the deterministic steady state is
+ model.getSteady(). It makes an approximation about it and stores the
+ derivatives to ‘rule_ders’ and ‘rule_ders_ss’. Also it runs a check() for
+ σ=0. */
+void
+ApproximationWelfare::approxAtSteady()
+{
+ welfare.calcDerivativesAtSteady();
+ 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
+
+ // construct the resulting decision rules
+ unc_fdr = std::make_unique<FoldDecisionRule>(*unc_ders, welfare.getModel().nys(), welfare.getModel().nexog(), welfare.getModel().getSteady());
+ 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
+{
+ KORD_RAISE_IF(!cond_fdr,
+ "Folded decision rule has not been created in ApproximationWelfare::getFoldCondWel");
+ return *cond_fdr;
+}
+
+void
+ApproximationWelfare::saveRuleDerivs(const FGSContainer &U, 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());
+ }
+}
diff --git a/mex/sources/k_order_welfare/approximation_welfare.hh b/mex/sources/k_order_welfare/approximation_welfare.hh
new file mode 100644
index 0000000000000000000000000000000000000000..c2d2d166f9e8c8d67933fac868772a367ed5a62b
--- /dev/null
+++ b/mex/sources/k_order_welfare/approximation_welfare.hh
@@ -0,0 +1,88 @@
+/*
+ * 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/>.
+ */
+
+#ifndef APPROXIMATION_WELFARE_H
+#define APPROXIMATION_WELFARE_H
+
+#include "k_ord_objective.hh"
+#include "journal.hh"
+
+#include <memory>
+
+class ApproximationWelfare
+{
+ KordwDynare &welfare;
+ 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);
+
+ const KordwDynare &
+ getWelfare() const
+ {
+ return welfare;
+ }
+ const FGSContainer &
+ get_rule_ders() const
+ {
+ return *rule_ders;
+ }
+ const FGSContainer &
+ get_rule_ders_s() const
+ {
+ 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();
+};
+
+#endif
diff --git a/mex/sources/k_order_welfare/k_ord_objective.cc b/mex/sources/k_order_welfare/k_ord_objective.cc
new file mode 100644
index 0000000000000000000000000000000000000000..cf5dbb67dbd02f8de920be93192b1fe796f39416
--- /dev/null
+++ b/mex/sources/k_order_welfare/k_ord_objective.cc
@@ -0,0 +1,389 @@
+/*
+ * 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/>.
+ */
+
+#include "k_ord_objective.hh"
+#include "objective_abstract_class.hh"
+
+#include <utility>
+#include <cassert>
+
+KordwDynare::KordwDynare(KordpDynare &m, ConstVector &NNZD_arg, Journal &jr, Vector &inParams, std::unique_ptr<ObjectiveAC> objectiveFile_arg, const std::vector<int> &dr_order) :
+ model{m}, NNZD{NNZD_arg}, journal{jr}, params{inParams}, resid(1), ud{1}, objectiveFile{std::move(objectiveFile_arg)}
+{
+ dynppToDyn = dr_order;
+ dynToDynpp.resize(model.ny());
+ for (int i = 0; i < model.ny(); i++)
+ dynToDynpp[dynppToDyn[i]] = i;
+}
+
+void
+KordwDynare::calcDerivativesAtSteady()
+{
+
+ assert(ud.begin() == ud.end());
+
+ std::vector<TwoDMatrix> dyn_ud; // Planner's objective derivatives, in Dynare form
+ dyn_ud.emplace_back(1, model.ny()); // Allocate Jacobian
+ dyn_ud.back().zeros();
+
+ for (int i = 2; i <= model.order(); i++)
+ {
+ // Higher order derivatives, as sparse (3-column) matrices
+ dyn_ud.emplace_back(static_cast<int>(NNZD[i-1]), 3);
+ dyn_ud.back().zeros();
+ }
+
+ Vector xx(model.nexog());
+ xx.zeros();
+ resid.zeros();
+ objectiveFile->eval(model.getSteady(), xx, params, resid, dyn_ud);
+
+ for (int i = 1; i <= model.order(); i++)
+ populateDerivativesContainer(dyn_ud, i);
+
+}
+
+void
+KordwDynare::populateDerivativesContainer(const std::vector<TwoDMatrix> &dyn_ud, int ord)
+{
+ const TwoDMatrix &u = dyn_ud[ord-1];
+
+ // utility derivatives FSSparseTensor instance
+ auto udTi = std::make_unique<FSSparseTensor>(ord, model.ny(), 1);
+
+ IntSequence s(ord, 0);
+
+ if (ord == 1)
+ for (int i = 0; i < u.ncols(); i++)
+ {
+ for (int j = 0; j < u.nrows(); j++)
+ {
+ double x = u.get(j, dynppToDyn[s[0]]);
+ if (x != 0.0)
+ udTi->insert(s, j, x);
+ }
+ s[0]++;
+ }
+ else // ord ≥ 2
+ for (int i = 0; i < u.nrows(); i++)
+ {
+ int j = static_cast<int>(u.get(i, 0))-1;
+ int i1 = static_cast<int>(u.get(i, 1))-1;
+ if (j < 0 || i1 < 0)
+ continue; // Discard empty entries (see comment in DynamicModelAC::unpackSparseMatrix())
+
+ for (int k = 0; k < ord; k++)
+ {
+ s[k] = dynToDynpp[i1 % model.ny()];
+ i1 /= model.ny();
+ }
+
+ if (ord == 2 && !s.isSorted())
+ continue; // Skip symmetric elements (only needed at order 2)
+ else if (ord > 2)
+ s.sort(); // For higher order, canonicalize the multi-index
+
+ double x = u.get(i, 2);
+ udTi->insert(s, j, x);
+ }
+
+ ud.insert(std::move(udTi));
+}
+
+template<>
+ctraits<Storage::unfold>::Tg &
+KOrderWelfare::g<Storage::unfold>()
+{
+ return _ug;
+}
+template<>
+const ctraits<Storage::unfold>::Tg &
+KOrderWelfare::g<Storage::unfold>() const
+{
+ return _ug;
+}
+template<>
+ctraits<Storage::fold>::Tg &
+KOrderWelfare::g<Storage::fold>()
+{
+ return _fg;
+}
+template<>
+const ctraits<Storage::fold>::Tg &
+KOrderWelfare::g<Storage::fold>() const
+{
+ return _fg;
+}
+template<>
+ctraits<Storage::unfold>::Tgs &
+KOrderWelfare::gs<Storage::unfold>()
+{
+ return _ugs;
+}
+template<>
+const ctraits<Storage::unfold>::Tgs &
+KOrderWelfare::gs<Storage::unfold>() const
+{
+ return _ugs;
+}
+template<>
+ctraits<Storage::fold>::Tgs &
+KOrderWelfare::gs<Storage::fold>()
+{
+ return _fgs;
+}
+template<>
+const ctraits<Storage::fold>::Tgs &
+KOrderWelfare::gs<Storage::fold>() const
+{
+ return _fgs;
+}
+
+template<>
+ctraits<Storage::unfold>::TU &
+KOrderWelfare::U<Storage::unfold>()
+{
+ return _uU;
+}
+template<>
+const ctraits<Storage::unfold>::TU &
+KOrderWelfare::U<Storage::unfold>() const
+{
+ return _uU;
+}
+template<>
+ctraits<Storage::fold>::TU &
+KOrderWelfare::U<Storage::fold>()
+{
+ return _fU;
+}
+template<>
+const ctraits<Storage::fold>::TU &
+KOrderWelfare::U<Storage::fold>() const
+{
+ return _fU;
+}
+template<>
+ctraits<Storage::unfold>::TW &
+KOrderWelfare::W<Storage::unfold>()
+{
+ return _uW;
+}
+template<>
+const ctraits<Storage::unfold>::TW &
+KOrderWelfare::W<Storage::unfold>() const
+{
+ return _uW;
+}
+template<>
+ctraits<Storage::fold>::TW &
+KOrderWelfare::W<Storage::fold>()
+{
+ return _fW;
+}
+template<>
+const ctraits<Storage::fold>::TW &
+KOrderWelfare::W<Storage::fold>() const
+{
+ return _fW;
+}
+template<>
+ctraits<Storage::unfold>::TWrond &
+KOrderWelfare::Wrond<Storage::unfold>()
+{
+ return _uWrond;
+}
+template<>
+const ctraits<Storage::unfold>::TWrond &
+KOrderWelfare::Wrond<Storage::unfold>() const
+{
+ return _uWrond;
+}
+template<>
+ctraits<Storage::fold>::TWrond &
+KOrderWelfare::Wrond<Storage::fold>()
+{
+ return _fWrond;
+}
+template<>
+const ctraits<Storage::fold>::TWrond &
+KOrderWelfare::Wrond<Storage::fold>() const
+{
+ return _fWrond;
+}
+template<>
+ctraits<Storage::unfold>::TGstack &
+KOrderWelfare::Gstack<Storage::unfold>()
+{
+ return _uGstack;
+}
+template<>
+const ctraits<Storage::unfold>::TGstack &
+KOrderWelfare::Gstack<Storage::unfold>() const
+{
+ return _uGstack;
+}
+template<>
+ctraits<Storage::fold>::TGstack &
+KOrderWelfare::Gstack<Storage::fold>()
+{
+ return _fGstack;
+}
+template<>
+const ctraits<Storage::fold>::TGstack &
+KOrderWelfare::Gstack<Storage::fold>() const
+{
+ return _fGstack;
+}
+template<>
+ctraits<Storage::unfold>::TXstack &
+KOrderWelfare::Xstack<Storage::unfold>()
+{
+ return _uXstack;
+}
+template<>
+const ctraits<Storage::unfold>::TXstack &
+KOrderWelfare::Xstack<Storage::unfold>() const
+{
+ return _uXstack;
+}
+template<>
+ctraits<Storage::fold>::TXstack &
+KOrderWelfare::Xstack<Storage::fold>()
+{
+ return _fXstack;
+}
+template<>
+const ctraits<Storage::fold>::TXstack &
+KOrderWelfare::Xstack<Storage::fold>() const
+{
+ return _fXstack;
+}
+template<>
+ctraits<Storage::unfold>::Tm &
+KOrderWelfare::m<Storage::unfold>()
+{
+ return _um;
+}
+template<>
+const ctraits<Storage::unfold>::Tm &
+KOrderWelfare::m<Storage::unfold>() const
+{
+ return _um;
+}
+template<>
+ctraits<Storage::fold>::Tm &
+KOrderWelfare::m<Storage::fold>()
+{
+ return _fm;
+}
+template<>
+const ctraits<Storage::fold>::Tm &
+KOrderWelfare::m<Storage::fold>() const
+{
+ return _fm;
+}
+
+KOrderWelfare::KOrderWelfare(int num_stat, int num_pred,
+ int num_both, int num_forw, int nu, int ord,
+ double discount_factor,
+ const TensorContainer<FSSparseTensor> &ucont,
+ const FGSContainer &g_arg,
+ const FGSContainer &gs_arg,
+ const TwoDMatrix &v,
+ Journal &jr) :
+ ypart(num_stat, num_pred, num_both, num_forw),
+ ny(ypart.ny()), nu(nu), maxk(ucont.getMaxDim()), order(ord),
+ discount_factor{discount_factor}, nvs{ypart.nys(), nu, nu, 1},
+ _uU(4), _fU(4), _uW(4), _fW(4), _uWrond(4), _fWrond(4), _ug(4),
+ _fg(g_arg), _ugs(4), _fgs(gs_arg), _uXstack(&_ug, ny),
+ _fXstack(&_fg, ny), _uGstack(&_ugs, ypart.nys(), nu),
+ _fGstack(&_fgs, ypart.nys(), nu), _um(maxk, v),
+ _fm(_um), u(ucont), journal(jr)
+{
+ KORD_RAISE_IF(v.ncols() != nu,
+ "Wrong number of columns of Vcov in KOrderWelfare constructor");
+ KORD_RAISE_IF(nu != v.nrows(),
+ "Wrong number of rows of Vcov in KOrderWelfare constructor");
+ for (int ord = 1; ord <= order; ord++)
+ {
+ JournalRecordPair pa(journal);
+ pa << u8"Unconditional welfare : performing step for order = " << ord << "\n" << endrec;
+ for (int j = 0; j <= ord; j++)
+ for (int i = 0; i <= j; i++)
+ {
+ Symmetry sym{ord-j, i, 0, j-i};
+ pa << "Recovering symmetry " << sym << "\n" << endrec;
+ auto U_sym = faaDiBrunoU<Storage::fold>(sym);
+ U<Storage::fold>().insert(std::move(U_sym));
+ }
+ }
+ U<Storage::unfold>() = UGSContainer(U<Storage::fold>());
+ g<Storage::unfold>() = UGSContainer(g<Storage::fold>());
+ gs<Storage::unfold>() = UGSContainer(gs<Storage::fold>());
+}
+
+// KOrderWelfare::sylvesterSolve() unfolded specialization
+/* Here we have an unfolded specialization of sylvesterSolve(). We simply
+ create the sylvester object and solve it. Note that the W_y is not
+ continuous in memory as assumed by the sylvester code, so we make a
+ temporary copy and pass it as matrix C.
+
+ If the B matrix is empty, in other words there are now forward looking
+ variables, then the system becomes AX=D which is solved by simple
+ matA.multInv().
+
+ If one wants to display the diagnostic messages from the Sylvester module,
+ then after the sylv.solve() one needs to call sylv.getParams().print(""). */
+
+template<>
+void
+KOrderWelfare::sylvesterSolve<Storage::unfold>(ctraits<Storage::unfold>::Ttensor &der) const
+{
+ JournalRecordPair pa(journal);
+ pa << "Sylvester equation for dimension = " << der.getSym()[0] << endrec;
+ KORD_RAISE_IF(!der.isFinite(),
+ "RHS of Sylverster is not finite");
+ TwoDMatrix gs_y(gs<Storage::unfold>().get(Symmetry{1, 0, 0, 0}));
+ TwoDMatrix A(1,1);
+ A.unit();
+ TwoDMatrix B(1,1);
+ B.unit();
+ B.mult(-discount_factor);
+ GeneralSylvester sylv(der.getSym()[0], 1, ypart.nys(),
+ 0,
+ A.getData(), B.getData(),
+ gs_y.getData(), der.getData());
+ sylv.solve();
+}
+
+// KOrder::sylvesterSolve() folded specialization
+/* Here is the folded specialization of sylvester. We unfold the right hand
+ side. Then we solve it by the unfolded version of sylvesterSolve(), and fold
+ it back and copy to output vector. */
+
+template<>
+void
+KOrderWelfare::sylvesterSolve<Storage::fold>(ctraits<Storage::fold>::Ttensor &der) const
+{
+ ctraits<Storage::unfold>::Ttensor tmp(der);
+ sylvesterSolve<Storage::unfold>(tmp);
+ ctraits<Storage::fold>::Ttensor ftmp(tmp);
+ der.getData() = const_cast<const Vector &>(ftmp.getData());
+}
diff --git a/mex/sources/k_order_welfare/k_ord_objective.hh b/mex/sources/k_order_welfare/k_ord_objective.hh
new file mode 100644
index 0000000000000000000000000000000000000000..66203bc7528247854b828ac81410b7f5308daec9
--- /dev/null
+++ b/mex/sources/k_order_welfare/k_ord_objective.hh
@@ -0,0 +1,653 @@
+/*
+ * 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/>.
+ */
+
+#ifndef K_ORD_OBJECTIVE_H
+#define K_ORD_OBJECTIVE_H
+
+#include "k_ord_dynare.hh"
+#include "objective_abstract_class.hh"
+
+class KordwDynare;
+
+class KordwDynare
+{
+public:
+ KordpDynare &model;
+ const ConstVector &NNZD;
+private:
+ Journal &journal;
+ Vector ¶ms;
+ Vector resid;
+ TensorContainer<FSSparseTensor> ud; // planner's objective derivatives, in Dynare++ form
+ std::vector<int> dynppToDyn; // Maps Dynare++ jacobian variable indices to Dynare ones
+ std::vector<int> dynToDynpp; // Maps Dynare jacobian variable indices to Dynare++ ones
+ std::unique_ptr<ObjectiveAC> objectiveFile;
+public:
+ KordwDynare(KordpDynare &m, ConstVector &NNZD_arg, Journal &jr, Vector &inParams, std::unique_ptr<ObjectiveAC> objectiveFile_arg, const std::vector<int> &varOrder);
+ void calcDerivativesAtSteady();
+ void populateDerivativesContainer(const std::vector<TwoDMatrix> &dyn_ud, int ord);
+ const TensorContainer<FSSparseTensor> &
+ getPlannerObjDerivatives() const
+ {
+ return ud;
+ }
+ const KordpDynare &
+ getModel() const
+ {
+ return model;
+ }
+ const Vector &
+ getResid() const
+ {
+ return resid;
+ }
+private:
+ /* Computes the permutations mapping back and forth between Dynare and
+ Dynare++ orderings of variables */
+ void computeJacobianPermutation(const std::vector<int> &var_order);
+
+};
+
+class KOrderWelfare
+{
+public:
+ const PartitionY ypart;
+ const int ny;
+ const int nu;
+ const int maxk;
+ const int order;
+ const double discount_factor;
+ /* Tensor variable dimension: variable sizes of all tensors in
+ containers, sizes of y, u and σ */
+ IntSequence nvs;
+ UGSContainer _uU;
+ FGSContainer _fU;
+ UGSContainer _uW;
+ FGSContainer _fW;
+ UGSContainer _uWrond;
+ FGSContainer _fWrond;
+ UGSContainer _ug;
+ FGSContainer _fg;
+ UGSContainer _ugs;
+ FGSContainer _fgs;
+ UnfoldedXContainer _uXstack;
+ FoldedXContainer _fXstack;
+ UnfoldedGContainer _uGstack;
+ FoldedGContainer _fGstack;
+
+ /* Moments: both folded and unfolded normal moment containers, both are
+ calculated at initialization */
+ UNormalMoments _um;
+ FNormalMoments _fm;
+
+ /* Planner objective derivatives: just a reference to the container of sparse
+ tensors of the derivatives, lives outside the class */
+ const TensorContainer<FSSparseTensor> &u;
+
+ /* These are the declarations of the template functions accessing the
+ containers. We declare template methods for accessing containers depending
+ on ‘fold’ and ‘unfold’ flag, we implement their specializations*/
+ template<Storage t>
+ typename ctraits<t>::Tg &g();
+ template<Storage t>
+ const typename ctraits<t>::Tg &g() const;
+ template<Storage t>
+ typename ctraits<t>::Tgs &gs();
+ template<Storage t>
+ const typename ctraits<t>::Tgs &gs() const;
+ template<Storage t>
+ typename ctraits<t>::TU &U();
+ template<Storage t>
+ const typename ctraits<t>::TU &U() const;
+ template<Storage t>
+ typename ctraits<t>::TW &W();
+ template<Storage t>
+ const typename ctraits<t>::TW &W() const;
+ template<Storage t>
+ typename ctraits<t>::TWrond &Wrond();
+ template<Storage t>
+ const typename ctraits<t>::TWrond &Wrond() const;
+ template<Storage t>
+ typename ctraits<t>::TXstack &Xstack();
+ template<Storage t>
+ const typename ctraits<t>::TXstack &Xstack() const;
+ template<Storage t>
+ typename ctraits<t>::TGstack &Gstack();
+ template<Storage t>
+ const typename ctraits<t>::TGstack &Gstack() const;
+ template<Storage t>
+ typename ctraits<t>::Tm &m();
+ template<Storage t>
+ const typename ctraits<t>::Tm &m() const;
+
+ Journal &journal;
+
+public:
+ KOrderWelfare(int num_stat, int num_pred, int num_both,
+ int num_forw, int nu, int ord, double discount_factor,
+ const TensorContainer<FSSparseTensor> &ucont,
+ const FGSContainer &g,
+ const FGSContainer &gs,
+ const TwoDMatrix &v,
+ Journal &jr);
+
+ /* Performs k-order step provided that k=2 or the k−1-th step has been
+ run, this is the core method */
+ template<Storage t>
+ void performStep(int order);
+ /* Calculates residuals of all solved equations for k-order and reports their
+ sizes, it is runnable after k-order performStep() has been run */
+ template<Storage t>
+ double check(int dim) const;
+
+ const FGSContainer &
+ getFoldU() const
+ {
+ return _fU;
+ }
+ const UGSContainer &
+ getUnfoldU() const
+ {
+ return _uU;
+ }
+ const FGSContainer &
+ getFoldW() const
+ {
+ return _fW;
+ }
+ const UGSContainer &
+ getUnfoldW() const
+ {
+ return _uW;
+ }
+ static bool
+ is_even(int i)
+ {
+ return i % 2 == 0;
+ }
+
+protected:
+ /* Calculates derivatives of U by Faà Di Bruno for the sparse container of
+ planner's objective derivatives and the container for the decision rule derivatives*/
+ template<Storage t>
+ std::unique_ptr<typename ctraits<t>::Ttensor> faaDiBrunoU(const Symmetry &sym) const;
+ /* Calculates derivatives of the compounded functions W and Gstack using the Faà Di Bruno formula*/
+ template<Storage t>
+ std::unique_ptr<typename ctraits<t>::Ttensor> faaDiBrunoW(const Symmetry &sym) const;
+
+ /* Solves the sylvester equation (templated fold, and unfold) */
+ template<Storage t>
+ void sylvesterSolve(typename ctraits<t>::Ttensor &der) const;
+
+ // Recovers W_y*ⁱ
+ template<Storage t>
+ void recover_y(int i);
+ // Recovers W_y*ⁱuʲ
+ template<Storage t>
+ void recover_yu(int i, int j);
+ // Recovers W_y*ⁱσʲ
+ template<Storage t>
+ void recover_ys(int i, int j);
+ // Recovers W_y*ⁱuʲσᵏ
+ template<Storage t>
+ void recover_yus(int i, int j, int k);
+ // Recovers W_σⁱ
+ template<Storage t>
+ void recover_s(int i);
+
+ // Calculates specified derivatives of Wrond and inserts them to the container
+ template<Storage t>
+ void fillWrond(int i, int j, int k);
+
+ // Calculates Dᵢⱼₖ
+ template<Storage t>
+ typename ctraits<t>::Ttensor calcH_ijk(int i, int j, int k) const;
+ template<Storage t>
+ typename ctraits<t>::Ttensor calcH_ik(int i, int k) const;
+ template<Storage t>
+ typename ctraits<t>::Ttensor calcH_k(int k) const;
+
+ // Calculates Eᵢⱼₖ
+ template<Storage t>
+ typename ctraits<t>::Ttensor calcJ_ijk(int i, int j, int k) const;
+ template<Storage t>
+ typename ctraits<t>::Ttensor calcJ_ik(int i, int k) const;
+ template<Storage t>
+ typename ctraits<t>::Ttensor calcJ_k(int k) const;
+
+};
+
+/* Here we implement Faà Di Bruno formula
+
+ ₖ ₗ
+ ∑ [u_zˡ]_γ₁…γₗ ∑ ∏ [h_{s^|cₘ|}]^γₘ
+ ˡ⁼¹ c∈ℳₗ,ₖ ᵐ⁼¹
+
+ where s is a given outer symmetry and k is the dimension of the symmetry. */
+template<Storage t>
+std::unique_ptr<typename ctraits<t>::Ttensor>
+KOrderWelfare::faaDiBrunoU(const Symmetry &sym) const
+{
+ JournalRecordPair pa(journal);
+ pa << u8"Faà Di Bruno U container for " << sym << endrec;
+ auto res = std::make_unique<typename ctraits<t>::Ttensor>(1, TensorDimens(sym, nvs));
+ FaaDiBruno bruno(journal);
+ bruno.calculate(Xstack<t>(), u, *res);
+ return res;
+}
+
+/* The same as KOrder::faaDiBrunoW(), but for W and G stack. */
+template<Storage t>
+std::unique_ptr<typename ctraits<t>::Ttensor>
+KOrderWelfare::faaDiBrunoW(const Symmetry &sym) const
+{
+ JournalRecordPair pa(journal);
+ pa << u8"Faà Di Bruno G container for " << sym << endrec;
+ TensorDimens tdims(sym, nvs);
+ auto res = std::make_unique<typename ctraits<t>::Ttensor>(1, tdims);
+ FaaDiBruno bruno(journal);
+ bruno.calculate(Gstack<t>(), W<t>(), *res);
+ return res;
+}
+
+template<Storage t>
+void
+KOrderWelfare::performStep(int order)
+{
+ KORD_RAISE_IF(order-1 != W<t>().getMaxDim() and order > 1, "Wrong order for KOrder::performStep");
+ JournalRecordPair pa(journal);
+ pa << "Performing step for order = " << order << endrec;
+
+ recover_y<t>(order);
+
+ for (int i = 0; i < order; i++)
+ recover_yu<t>(i, order-i);
+
+ recover_ys<t>(order-1, 1);
+
+ for (int j = 2; j < order; j++)
+ {
+ for (int i = 1; i <= j-1; i++)
+ recover_yus<t>(order-j, i, j-i);
+ recover_ys<t>(order-j, j);
+ recover_yus<t>(0, order-j, j);
+ }
+
+ recover_s<t>(order);
+}
+
+/* Here we solve [F_yⁱ]=0. First we calculate conditional W_yⁱ (it misses l=i since W_yⁱ does not exist yet). Then calculate conditional F_yⁱ and
+ we have the right hand side of equation. Since we miss two orders, we solve
+ by Sylvester, and the solution as the derivative g_yⁱ. Then we need
+ to update G_yⁱ running multAndAdd() for both dimensions 1 and i.
+
+ Requires: everything at order ≤ i−1
+
+ Provides: W_yⁱ and Wrond_yⁱ
+*/
+template<Storage t>
+void
+KOrderWelfare::recover_y(int i)
+{
+ Symmetry sym{i, 0, 0, 0};
+ JournalRecordPair pa(journal);
+ pa << "Conditional welfare: recovering symmetry " << sym << "\n" << endrec;
+
+ auto Wrond_yi = faaDiBrunoW<t>(sym);
+ auto Wrond_yi_ptr = Wrond_yi.get();
+ Wrond<t>().insert(std::move(Wrond_yi));
+
+ auto W_yi = U<t>().get(sym);
+ W_yi.add(discount_factor, *Wrond_yi_ptr);
+
+ sylvesterSolve<t>(W_yi);
+
+ W<t>().insert(std::make_unique<typename ctraits<t>::Ttensor>(W_yi));
+
+ Gstack<t>().multAndAdd(i, W<t>(), *Wrond_yi_ptr);
+}
+
+/* Here we solve [F_yⁱuʲ]=0 to obtain W_yⁱuʲ for j>0. We calculate conditional
+ Wrond_yⁱuʲ and calculate conditional F_yⁱuʲ and we have
+ the right hand side.
+
+ Requires: everything at order ≤ i+j−1, Wrond_yⁱ⁺ʲ and W_yⁱ⁺ʲ.
+
+ Provides: W_yⁱuʲ and Wrond_yⁱuʲ
+*/
+template<Storage t>
+void
+KOrderWelfare::recover_yu(int i, int j)
+{
+ Symmetry sym{i, j, 0, 0};
+ JournalRecordPair pa(journal);
+ pa << "Conditional welfare: recovering symmetry " << sym << endrec;
+
+ auto Wrond_yiuj = faaDiBrunoW<t>(sym);
+ auto Wrond_yiuj_ptr = Wrond_yiuj.get();
+ Wrond<t>().insert(std::move(Wrond_yiuj));
+
+ auto W_yiuj = U<t>().get(sym);
+ W_yiuj.add(discount_factor, *Wrond_yiuj_ptr);
+
+ W<t>().insert(std::make_unique<typename ctraits<t>::Ttensor>(W_yiuj));
+
+}
+
+/* Here we solve [F_yⁱσʲ]+[Dᵢⱼ]+[Eᵢⱼ]=0 to obtain W_yⁱσʲ. We calculate
+ conditional Wrond_yⁱσʲ (missing dimensions i+j), calculate conditional
+ F_yⁱσʲ. Before we can calculate Dᵢⱼ and Eᵢⱼ, we have to calculate
+ Wrond_yⁱu′ᵐσʲ⁻ᵐ for m=1,…,j. Then we add the Dᵢⱼ and Eᵢⱼ to obtain the right
+ hand side. Then we solve the sylvester to obtain g_yⁱσʲ. Then we update
+ G_yⁱσʲ for l=1 and l=i+j.
+
+ Requires: everything at order ≤ i+j−1, g_yⁱ⁺ʲ, G_yⁱu′ʲ and g_yⁱuʲ through
+ Dᵢⱼ, g_yⁱuᵐσʲ⁻ᵐ for m=1,…,j−1 through Eᵢⱼ.
+
+ Provides: g_yⁱσʲ and G_yⁱσʲ, and finally G_yⁱu′ᵐσʲ⁻ᵐ for m=1,…,j. The latter
+ is calculated by fillG() before the actual calculation.
+*/
+template<Storage t>
+void
+KOrderWelfare::recover_ys(int i, int j)
+{
+ Symmetry sym{i, 0, 0, j};
+ JournalRecordPair pa(journal);
+ pa << "Conditional welfare: recovering symmetry " << sym << endrec;
+
+ fillWrond<t>(i, 0, j);
+
+ if (is_even(j))
+ {
+ auto Wrond_yisj = faaDiBrunoW<t>(sym);
+ auto Wrond_yisj_ptr = Wrond_yisj.get();
+ Wrond<t>().insert(std::move(Wrond_yisj));
+
+ auto W_yisj = U<t>().get(sym);
+ W_yisj.add(discount_factor, *Wrond_yisj_ptr);
+ {
+ auto H_ij = calcH_ik<t>(i, j);
+ W_yisj.add(-1.0, H_ij);
+ }
+
+ if (j >= 3)
+ {
+ auto J_ij = calcJ_ik<t>(i, j);
+ W_yisj.add(-1.0, J_ij);
+ }
+
+ sylvesterSolve<t>(W_yisj);
+
+ W<t>().insert(std::make_unique<typename ctraits<t>::Ttensor>(W_yisj));
+
+ Gstack<t>().multAndAdd(i+j, W<t>(), *Wrond_yisj_ptr);
+ }
+}
+
+/* Here we solve [F_yⁱuʲσᵏ]+[Dᵢⱼₖ]+[Eᵢⱼₖ]=0 to obtain g_yⁱuʲσᵏ. First we
+ calculate conditional G_yⁱuʲσᵏ (missing only for dimension l=1), then we
+ evaluate conditional F_yⁱuʲσᵏ. Before we can calculate Dᵢⱼₖ, and Eᵢⱼₖ, we
+ need to G_yⁱuʲu′ᵐσᵏ⁻ᵐ for m=1,…,k. This is done by fillG(). Then we
+ have right hand side and we multiply by A⁻¹ to obtain g_yⁱuʲσᵏ. Finally we
+ have to update G_yⁱuʲσᵏ by multAndAdd() for dimension l=1.
+
+ Requires: everything at order ≤ i+j+k, g_yⁱ⁺ʲσᵏ through G_yⁱuʲσᵏ involved in
+ right hand side, then g_yⁱuʲ⁺ᵏ through Dᵢⱼₖ, and g_yⁱuʲ⁺ᵐσᵏ⁻ᵐ for m=1,…,k−1
+ through Eᵢⱼₖ.
+
+ Provides: g_yⁱuʲσᵏ, G_yⁱuʲσᵏ, and G_yⁱuʲu′ᵐσᵏ⁻ᵐ for m=1,…,k
+*/
+template<Storage t>
+void
+KOrderWelfare::recover_yus(int i, int j, int k)
+{
+ Symmetry sym{i, j, 0, k};
+ JournalRecordPair pa(journal);
+ pa << "Conditional welfare: recovering symmetry " << sym << endrec;
+
+ fillWrond<t>(i, j, k);
+
+ if (is_even(k))
+ {
+ auto Wrond_yiujsk = faaDiBrunoW<t>(sym);
+ auto Wrond_yiujsk_ptr = Wrond_yiujsk.get();
+ Wrond<t>().insert(std::move(Wrond_yiujsk));
+
+ auto W_yiujsk = U<t>().get(sym);
+ W_yiujsk.add(discount_factor, *Wrond_yiujsk_ptr);
+
+ auto H_ijk = calcH_ijk<t>(i, j, k);
+ W_yiujsk.add(-1.0, H_ijk);
+
+ if (k >= 3)
+ {
+ auto J_ijk = calcJ_ijk<t>(i, j, k);
+ W_yiujsk.add(-1.0, J_ijk);
+ }
+
+ W<t>().insert(std::make_unique<typename ctraits<t>::Ttensor>(W_yiujsk));
+
+ Gstack<t>().multAndAdd(i+j+k, W<t>(), *Wrond_yiujsk_ptr);
+ }
+}
+
+/* Here we solve [F_{σⁱ}]+[Dᵢ]+[Eᵢ]=0 to recover g_σⁱ. First we calculate
+ conditional Wrond_σⁱ (missing dimension and l=i), then we calculate
+ conditional F_σⁱ. Before we can calculate Dᵢ and Eᵢ, we have to obtain
+ G_u′ᵐσⁱ⁻ᵐ for m=1,…,i. Then adding Dᵢ and Eᵢ we have the right hand side. We
+ solve by S⁻¹ multiplication and update G_σⁱ by calling multAndAdd() for
+ dimension l=1.
+
+ Recall that the solved equation here is:
+ ny(ny), nu(nu), F_u′ʲσⁱ⁻ʲ.
+
+ Provides: g_σⁱ, G_σⁱ, and G_u′ᵐσⁱ⁻ᵐ for m=1,…,i
+*/
+template<Storage t>
+void
+KOrderWelfare::recover_s(int i)
+{
+ Symmetry sym{0, 0, 0, i};
+ JournalRecordPair pa(journal);
+ pa << "Conditional welfare: recovering symmetry " << sym << endrec;
+
+ fillWrond<t>(0, 0, i);
+
+ if (is_even(i))
+ {
+ auto Wrond_si = faaDiBrunoW<t>(sym);
+ auto Wrond_si_ptr = Wrond_si.get();
+ Wrond<t>().insert(std::move(Wrond_si));
+
+ auto W_si = U<t>().get(sym);
+ W_si.add(discount_factor, *Wrond_si_ptr);
+
+ {
+ auto H_i = calcH_k<t>(i);
+ W_si.add(-1.0, H_i);
+ }
+
+ if (i >= 3)
+ {
+ auto J_i = calcJ_k<t>(i);
+ W_si.add(-1.0, J_i);
+ }
+
+ W_si.mult(1/(1-discount_factor));
+
+ W<t>().insert(std::make_unique<typename ctraits<t>::Ttensor>(W_si));
+
+ Gstack<t>().multAndAdd(i, W<t>(), *Wrond_si_ptr);
+ }
+}
+
+/* Here we calculate and insert Wrond_yⁱuʲu′ᵐσᵏ⁻ᵐ for m=1,…,k. The derivatives are inserted only for k−m being even. */
+template<Storage t>
+void
+KOrderWelfare::fillWrond(int i, int j, int k)
+{
+ for (int m = 1; m <= k; m++)
+ if (is_even(k-m))
+ {
+ auto Wrond_yiujupms = faaDiBrunoW<t>(Symmetry{i, j, m, k-m});
+ Wrond<t>().insert(std::move(Wrond_yiujupms));
+ }
+}
+
+/* Here we calculate:
+
+ [Hᵢⱼₖ]_α₁…αᵢβ₁…βⱼ = [F_yⁱuʲu′ᵏ]_α₁…αᵢβ₁…βⱼγ₁…γₖ [Σ]^γ₁…γₖ
+ = -β [Wrond_yⁱuʲu′ᵏ]_α₁…αᵢβ₁…βⱼγ₁…γₖ [Σ]^γ₁…γₖ
+ So it is non zero only for even k. */
+template<Storage t>
+typename ctraits<t>::Ttensor
+KOrderWelfare::calcH_ijk(int i, int j, int k) const
+{
+ typename ctraits<t>::Ttensor res(1, TensorDimens(Symmetry{i, j, 0, 0}, nvs));
+ res.zeros();
+ if (is_even(k))
+ {
+ auto tmp = faaDiBrunoW<t>(Symmetry{i, j, k, 0});
+ tmp->contractAndAdd(2, res, m<t>().get(Symmetry{k}));
+ }
+ res.mult(-discount_factor);
+ return res;
+}
+
+/* Here we calculate
+ ₖ₋₁ ⎛k⎞
+ [Jᵢⱼₖ]_α₁…αᵢβ₁…βⱼ = ∑ ⎝m⎠ [F_yⁱuʲu′ᵐσᵏ⁻ᵐ]_α₁…αᵢβ₁…βⱼγ₁…γₘ [Σ]^γ₁…γₘ
+ ᵐ⁼¹
+ ₖ₋₁ ⎛k⎞
+ = -β ∑ ⎝m⎠ [Wrond_yⁱuʲu′ᵐσᵏ⁻ᵐ]_α₁…αᵢβ₁…βⱼγ₁…γₘ [Σ]^γ₁…γₘ
+ ᵐ⁼¹
+The sum can sum only for even m. */
+template<Storage t>
+typename ctraits<t>::Ttensor
+KOrderWelfare::calcJ_ijk(int i, int j, int k) const
+{
+ typename ctraits<t>::Ttensor res(1, TensorDimens(Symmetry{i, j, 0, 0}, nvs));
+ res.zeros();
+ for (int n = 2; n <= k-1; n += 2)
+ {
+ auto tmp = faaDiBrunoW<t>(Symmetry{i, j, n, k-n});
+ tmp->mult(static_cast<double>(PascalTriangle::noverk(k, n)));
+ tmp->contractAndAdd(2, res, m<t>().get(Symmetry{n}));
+ }
+ res.mult(-discount_factor);
+ return res;
+}
+
+template<Storage t>
+typename ctraits<t>::Ttensor
+KOrderWelfare::calcH_ik(int i, int k) const
+{
+ return calcH_ijk<t>(i, 0, k);
+}
+
+template<Storage t>
+typename ctraits<t>::Ttensor
+KOrderWelfare::calcH_k(int k) const
+{
+ return calcH_ijk<t>(0, 0, k);
+}
+
+template<Storage t>
+typename ctraits<t>::Ttensor
+KOrderWelfare::calcJ_ik(int i, int k) const
+{
+ return calcJ_ijk<t>(i, 0, k);
+}
+
+template<Storage t>
+typename ctraits<t>::Ttensor
+KOrderWelfare::calcJ_k(int k) const
+{
+ return calcJ_ijk<t>(0, 0, k);
+}
+
+/* Here we check for residuals of all the solved equations at the given order.
+ The method returns the largest residual size. Each check simply evaluates
+ the equation. */
+
+template<Storage t>
+double
+KOrderWelfare::check(int dim) const
+{
+ KORD_RAISE_IF(dim > W<t>().getMaxDim(),
+ "Wrong dimension for KOrderWelfare::check");
+ JournalRecordPair pa(journal);
+ pa << "Checking residuals for order = " << dim << endrec;
+
+ double maxerror = 0.0;
+
+ // Check for F_yⁱuʲ=0
+ for (int i = 0; i <= dim; i++)
+ {
+ Symmetry sym{dim-i, i, 0, 0};
+ auto r = W<t>().get(sym);
+ r.add(-1.0, U<t>().get(sym));
+ r.add(-discount_factor, Wrond<t>().get(sym));
+ double err = r.getData().getMax();
+ JournalRecord(journal) << "\terror for symmetry " << sym << "\tis " << err << endrec;
+ maxerror = std::max(err, maxerror);
+ }
+
+ // Check for F_yⁱuʲu′ᵏ+Hᵢⱼₖ+Jᵢⱼₖ=0
+ for (auto &si : SymmetrySet(dim, 3))
+ {
+ int i = si[0];
+ int j = si[1];
+ int k = si[2];
+ if (i+j > 0 && k > 0 && is_even(k))
+ {
+ Symmetry sym{i, j, 0, k};
+ auto r = W<t>().get(sym);
+ r.add(-1.0, U<t>().get(sym));
+ r.add(-discount_factor, Wrond<t>().get(sym));
+ auto H_ijk = calcH_ijk<t>(i, j, k);
+ r.add(1.0, H_ijk);
+ auto J_ijk = calcJ_ijk<t>(i, j, k);
+ r.add(1.0, J_ijk);
+ double err = r.getData().getMax();
+ JournalRecord(journal) << "\terror for symmetry " << sym << "\tis " << err << endrec;
+ maxerror = std::max(err, maxerror);
+ }
+ }
+
+ // Check for F_σⁱ+Dᵢ+Eᵢ=0
+ if (is_even(dim))
+ {
+ Symmetry sym{0, 0, 0, dim};
+ auto r = W<t>().get(sym);
+ r.add(-1.0, U<t>().get(sym));
+ r.add(-discount_factor, Wrond<t>().get(sym));
+ auto H_k = calcH_k<t>(dim);
+ r.add(1.0, H_k);
+ auto J_k = calcJ_k<t>(dim);
+ r.add(1.0, J_k);
+ double err = r.getData().getMax();
+ JournalRecord(journal) << "\terror for symmetry " << sym << "\tis " << err << endrec;
+ maxerror = std::max(err, maxerror);
+ }
+
+ return maxerror;
+}
+
+#endif
diff --git a/mex/sources/k_order_welfare/k_order_welfare.cc b/mex/sources/k_order_welfare/k_order_welfare.cc
new file mode 100644
index 0000000000000000000000000000000000000000..8a41b91d3cb89ea8117346c7bdb881dd92ac5bce
--- /dev/null
+++ b/mex/sources/k_order_welfare/k_order_welfare.cc
@@ -0,0 +1,392 @@
+/*
+ * 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/>.
+ */
+
+#include "dynamic_m.hh"
+#include "dynamic_dll.hh"
+#include "objective_m.hh"
+
+#include "approximation.hh"
+#include "approximation_welfare.hh"
+#include "exception.hh"
+#include "dynare_exception.hh"
+#include "kord_exception.hh"
+#include "tl_exception.hh"
+#include "SylvException.hh"
+
+#include <algorithm>
+#include <cassert>
+
+#include "dynmex.h"
+
+/* Convert MATLAB Dynare endo and exo names cell array to a vector<string> array of
+ string pointers. */
+std::vector<std::string>
+DynareMxArrayToString(const mxArray *mxFldp)
+{
+ assert(mxIsCell(mxFldp));
+ std::vector<std::string> r;
+ for (size_t i = 0; i < mxGetNumberOfElements(mxFldp); i++)
+ {
+ const mxArray *cell_mx = mxGetCell(mxFldp, i);
+ if (!(cell_mx && mxIsChar(cell_mx)))
+ mexErrMsgTxt("Cell is not a character array");
+ r.emplace_back(mxArrayToString(cell_mx));
+ }
+
+ return r;
+}
+
+/* Vector for storing field names like “W_0”, “W_1”, …
+ A static structure is needed since MATLAB apparently does not create its own
+ copy of the strings (contrary to what is said at:
+ https://fr.mathworks.com/matlabcentral/answers/315937-mxcreatestructarray-and-mxcreatestructmatrix-field-name-memory-management
+ )
+*/
+
+std::vector<std::string> U_fieldnames;
+std::vector<std::string> W_fieldnames;
+
+void
+copy_derivatives(mxArray *destin, const Symmetry &sym, const FGSContainer &derivs, const char *fieldname)
+{
+ const FGSTensor &x = derivs.get(sym);
+ auto x_unfolded = x.unfold();
+ int n = x_unfolded->nrows();
+ int m = x_unfolded->ncols();
+ mxArray *tmp = mxCreateDoubleMatrix(n, m, mxREAL);
+ std::copy_n(x_unfolded->getData().base(), n*m, mxGetPr(tmp));
+ mxSetField(destin, 0, fieldname, tmp);
+}
+
+/* The k_order_welfare function computes the conditional welfare starting from the non-stochastic steady state and the derivatives of the felicity and welfare functions U(h(yₜ₋₁, uₜ, σ)) and W(yₜ₋₁, uₜ, σ). The unconditional welfare shall be introduced in a future version with the update of the dynare_simul_.
+
+The routine proceeds in several steps:
+1. Computation of the kth-order policy functions: the code is a copy-paste from the k_order_perturbation.cc file
+2. Computation of the kth-order derivatives of the felicity and of the welfare functions. The call to the approxAtSteady method in the ApproximationWelfare class carries out the necessary operation
+ - Importing the derivatives of the felicity function with the calcDerivativesAtSteady() method of the KordwDynare class. It relies on the Matlab-generated files, which are handled by the ObjectiveAC and ObjectiveMFile classes
+ - Pinpointing the derivatives of the felicity and welfare functions. The performStep method of the KOrderWelfare class carries out the calculations,resorting to the FaaDiBruno class and its methods to get the needed intermediary results.
+*/
+
+extern "C" {
+
+ void mexFunction(int nlhs, mxArray *plhs[],
+ int nrhs, const mxArray *prhs[])
+ {
+
+ mexPrintf("k_order_welfare\n");
+
+ const mxArray *dr_mx = prhs[0];
+ const mxArray *M_mx = prhs[1];
+ const mxArray *options_mx = prhs[2];
+
+ auto get_int_field = [](const mxArray *struct_mx, const std::string &fieldname)
+ {
+ mxArray *field_mx = mxGetField(struct_mx, 0, fieldname.c_str());
+ if (!(field_mx && mxIsScalar(field_mx) && mxIsNumeric(field_mx)))
+ mexErrMsgTxt(("Field `" + fieldname + "' should be a numeric scalar").c_str());
+ return static_cast<int>(mxGetScalar(field_mx));
+ };
+
+ // Extract various fields from options_
+ const mxArray *ramsey_policy_mx = mxGetField(options_mx, 0, "ramsey_policy");
+ if (!(ramsey_policy_mx && mxIsLogicalScalar(ramsey_policy_mx)))
+ mexErrMsgTxt("options_.ramsey_policy should be a logical scalar");
+ bool ramsey_policy = static_cast<bool>(mxGetScalar(ramsey_policy_mx));
+
+ if (ramsey_policy == false)
+ mexErrMsgTxt("The considered model must be a Ramsey-typed model !");
+
+
+ const int kOrder = get_int_field(options_mx, "order");
+ if (kOrder < 1)
+ mexErrMsgTxt("options_.order must be at least 1");
+
+ const mxArray *use_dll_mx = mxGetField(options_mx, 0, "use_dll");
+ if (!(use_dll_mx && mxIsLogicalScalar(use_dll_mx)))
+ mexErrMsgTxt("options_.use_dll should be a logical scalar");
+ bool use_dll = static_cast<bool>(mxGetScalar(use_dll_mx));
+
+ double qz_criterium = 1+1e-6;
+ const mxArray *qz_criterium_mx = mxGetField(options_mx, 0, "qz_criterium");
+ if (qz_criterium_mx && mxIsScalar(qz_criterium_mx) && mxIsNumeric(qz_criterium_mx))
+ qz_criterium = mxGetScalar(qz_criterium_mx);
+
+ const mxArray *threads_mx = mxGetField(options_mx, 0, "threads");
+ if (!threads_mx)
+ mexErrMsgTxt("Can't find field options_.threads");
+ const mxArray *num_threads_mx = mxGetField(threads_mx, 0, "k_order_perturbation");
+ if (!(num_threads_mx && mxIsScalar(num_threads_mx) && mxIsNumeric(num_threads_mx)))
+ mexErrMsgTxt("options_.threads.k_order_perturbation be a numeric scalar");
+ int num_threads = static_cast<int>(mxGetScalar(num_threads_mx));
+
+ const mxArray *debug_mx = mxGetField(options_mx, 0, "debug");
+ if (!(debug_mx && mxIsLogicalScalar(debug_mx)))
+ mexErrMsgTxt("options_.debug should be a logical scalar");
+ bool debug = static_cast<bool>(mxGetScalar(debug_mx));
+
+ // Extract various fields from M_
+ const mxArray *fname_mx = mxGetField(M_mx, 0, "fname");
+ if (!(fname_mx && mxIsChar(fname_mx) && mxGetM(fname_mx) == 1))
+ mexErrMsgTxt("M_.fname should be a character string");
+ std::string fname{mxArrayToString(fname_mx)};
+
+ const mxArray *params_mx = mxGetField(M_mx, 0, "params");
+ if (!(params_mx && mxIsDouble(params_mx)))
+ mexErrMsgTxt("M_.params should be a double precision array");
+ Vector modParams{ConstVector{params_mx}};
+ if (!modParams.isFinite())
+ mexErrMsgTxt("M_.params contains NaN or Inf");
+
+ const mxArray *param_names_mx = mxGetField(M_mx, 0, "param_names");
+ if (!(param_names_mx && mxIsCell(param_names_mx)))
+ mexErrMsgTxt("M_.param_names should be a cell array");
+ std::vector<std::string> paramNames = DynareMxArrayToString(param_names_mx);
+ auto it = std::find(paramNames.begin(), paramNames.end(), "optimal_policy_discount_factor");
+ double discount_factor;
+ if (it != paramNames.end())
+ discount_factor = modParams[std::distance(paramNames.begin(), it)];
+ else
+ {
+ mexErrMsgTxt("M_.param_names does not contain any \"optimal_policy_discount_factor\".");
+ return; // To silence a GCC warning about discount_factor unitialized
+ }
+
+ const mxArray *sigma_e_mx = mxGetField(M_mx, 0, "Sigma_e");
+ if (!(sigma_e_mx && mxIsDouble(sigma_e_mx) && mxGetM(sigma_e_mx) == mxGetN(sigma_e_mx)))
+ mexErrMsgTxt("M_.Sigma_e should be a double precision square matrix");
+ TwoDMatrix vCov{ConstTwoDMatrix{sigma_e_mx}};
+ if (!vCov.isFinite())
+ mexErrMsgTxt("M_.Sigma_e contains NaN or Inf");
+
+ const int nStat = get_int_field(M_mx, "nstatic");
+ const int nPred = get_int_field(M_mx, "npred");
+ const int nBoth = get_int_field(M_mx, "nboth");
+ const int nForw = get_int_field(M_mx, "nfwrd");
+
+ const int nExog = get_int_field(M_mx, "exo_nbr");
+ const int nEndo = get_int_field(M_mx, "endo_nbr");
+ const int nPar = get_int_field(M_mx, "param_nbr");
+
+ const mxArray *lead_lag_incidence_mx = mxGetField(M_mx, 0, "lead_lag_incidence");
+ if (!(lead_lag_incidence_mx && mxIsDouble(lead_lag_incidence_mx) && mxGetM(lead_lag_incidence_mx) == 3
+ && mxGetN(lead_lag_incidence_mx) == static_cast<size_t>(nEndo)))
+ mexErrMsgTxt("M_.lead_lag_incidence should be a double precision matrix with 3 rows and M_.endo_nbr columns");
+ ConstTwoDMatrix llincidence{lead_lag_incidence_mx};
+
+ const mxArray *nnzderivatives_mx = mxGetField(M_mx, 0, "NNZDerivatives");
+ if (!(nnzderivatives_mx && mxIsDouble(nnzderivatives_mx)))
+ mexErrMsgTxt("M_.NNZDerivatives should be a double precision array");
+ ConstVector NNZD{nnzderivatives_mx};
+ if (NNZD.length() < kOrder || NNZD[kOrder-1] == -1)
+ mexErrMsgTxt("The derivatives were not computed for the required order. Make sure that you used the right order option inside the `stoch_simul' command");
+
+ const mxArray *nnzderivatives_obj_mx = mxGetField(M_mx, 0, "NNZDerivatives_objective");
+ if (!(nnzderivatives_obj_mx && mxIsDouble(nnzderivatives_obj_mx)))
+ mexErrMsgTxt("M_.NNZDerivatives should be a double precision array");
+ ConstVector NNZD_obj{nnzderivatives_obj_mx};
+ if (NNZD.length() < kOrder || NNZD_obj[kOrder-1] == -1)
+ mexErrMsgTxt("The derivatives were not computed for the required order. Make sure that you used the right order option inside the `stoch_simul' command");
+
+ const mxArray *endo_names_mx = mxGetField(M_mx, 0, "endo_names");
+ if (!(endo_names_mx && mxIsCell(endo_names_mx) && mxGetNumberOfElements(endo_names_mx) == static_cast<size_t>(nEndo)))
+ mexErrMsgTxt("M_.endo_names should be a cell array of M_.endo_nbr elements");
+ std::vector<std::string> endoNames = DynareMxArrayToString(endo_names_mx);
+
+ const mxArray *exo_names_mx = mxGetField(M_mx, 0, "exo_names");
+ if (!(exo_names_mx && mxIsCell(exo_names_mx) && mxGetNumberOfElements(exo_names_mx) == static_cast<size_t>(nExog)))
+ mexErrMsgTxt("M_.exo_names should be a cell array of M_.exo_nbr elements");
+ std::vector<std::string> exoNames = DynareMxArrayToString(exo_names_mx);
+
+ const mxArray *dynamic_tmp_nbr_mx = mxGetField(M_mx, 0, "dynamic_tmp_nbr");
+ if (!(dynamic_tmp_nbr_mx && mxIsDouble(dynamic_tmp_nbr_mx) && mxGetNumberOfElements(dynamic_tmp_nbr_mx) >= static_cast<size_t>(kOrder+1)))
+ mexErrMsgTxt("M_.dynamic_tmp_nbr should be a double precision array with strictly more elements than the order of derivation");
+ int ntt = std::accumulate(mxGetPr(dynamic_tmp_nbr_mx), mxGetPr(dynamic_tmp_nbr_mx)+kOrder+1, 0);
+
+ // Extract various fields from dr
+ const mxArray *ys_mx = mxGetField(dr_mx, 0, "ys"); // and not in order of dr.order_var
+ if (!(ys_mx && mxIsDouble(ys_mx)))
+ mexErrMsgTxt("dr.ys should be a double precision array");
+ Vector ySteady{ConstVector{ys_mx}};
+ if (!ySteady.isFinite())
+ mexErrMsgTxt("dr.ys contains NaN or Inf");
+
+ const mxArray *order_var_mx = mxGetField(dr_mx, 0, "order_var");
+ if (!(order_var_mx && mxIsDouble(order_var_mx) && mxGetNumberOfElements(order_var_mx) == static_cast<size_t>(nEndo)))
+ mexErrMsgTxt("dr.order_var should be a double precision array of M_.endo_nbr elements");
+ std::vector<int> dr_order(nEndo);
+ std::transform(mxGetPr(order_var_mx), mxGetPr(order_var_mx)+nEndo, dr_order.begin(),
+ [](double x) { return static_cast<int>(x)-1; });
+
+ const int nSteps = 0; // Dynare++ solving steps, for time being default to 0 = deterministic steady state
+
+ // Journal is not written on-disk, unless options_.debug = true (see #1735)
+ Journal journal;
+ if (debug)
+ journal = Journal{fname + ".jnl"};
+
+ std::unique_ptr<DynamicModelAC> dynamicModelFile;
+ if (use_dll)
+ dynamicModelFile = std::make_unique<DynamicModelDLL>(fname, ntt, kOrder);
+ else
+ dynamicModelFile = std::make_unique<DynamicModelMFile>(fname, ntt);
+
+ // intiate tensor library
+ TLStatic::init(kOrder, nStat+2*nPred+3*nBoth+2*nForw+nExog);
+
+ // Set number of parallel threads
+ sthread::detach_thread_group::max_parallel_threads = num_threads;
+
+ // make KordpDynare object
+ KordpDynare dynare(endoNames, exoNames, nExog, nPar,
+ ySteady, vCov, modParams, nStat, nPred, nForw, nBoth,
+ NNZD, nSteps, kOrder, journal, std::move(dynamicModelFile),
+ dr_order, llincidence);
+
+
+ // construct main K-order approximation class
+ Approximation app(dynare, journal, nSteps, false, qz_criterium);
+ // run stochastic steady
+ app.walkStochSteady();
+
+ const mxArray *objective_tmp_nbr_mx = mxGetField(M_mx, 0, "objective_tmp_nbr");
+ if (!(objective_tmp_nbr_mx && mxIsDouble(objective_tmp_nbr_mx) && mxGetNumberOfElements(objective_tmp_nbr_mx) >= static_cast<size_t>(kOrder+1)))
+ mexErrMsgTxt("M_.objective_tmp_nbr should be a double precision array with strictly more elements than the order of derivation");
+ int ntt_objective = std::accumulate(mxGetPr(objective_tmp_nbr_mx), mxGetPr(objective_tmp_nbr_mx)+kOrder+1, 0);
+
+ //Getting derivatives of the planner's objective function
+ std::unique_ptr<ObjectiveAC> objectiveFile;
+ objectiveFile = std::make_unique<ObjectiveMFile>(fname, ntt_objective);
+
+ // make KordwDynare object
+ KordwDynare welfare(dynare, NNZD_obj, journal, modParams, std::move(objectiveFile), dr_order);
+
+ // construct main K-order approximation class of welfare
+ 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");
+ }
+ }
+ }
+ }
+ }
+ } // end of mexFunction()
+} // end of extern C
diff --git a/mex/sources/k_order_welfare/objective_abstract_class.hh b/mex/sources/k_order_welfare/objective_abstract_class.hh
new file mode 100644
index 0000000000000000000000000000000000000000..868cf129fab1befe671c9c9fcaf54993dc5e158c
--- /dev/null
+++ b/mex/sources/k_order_welfare/objective_abstract_class.hh
@@ -0,0 +1,38 @@
+/*
+ * 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/>.
+ */
+
+#ifndef _OBJECTIVEAC_HH
+#define _OBJECTIVEAC_HH
+
+#include <vector>
+
+#include "twod_matrix.hh"
+
+class ObjectiveAC
+{
+protected:
+ int ntt; // Size of vector of temporary terms
+public:
+ ObjectiveAC(int ntt_arg) : ntt{ntt_arg}
+ {
+ };
+ virtual ~ObjectiveAC() = default;
+ virtual void eval(const Vector &y, const Vector &x, const Vector ¶ms, Vector &residual, std::vector<TwoDMatrix> &md) = 0;
+};
+#endif
diff --git a/mex/sources/k_order_welfare/objective_m.cc b/mex/sources/k_order_welfare/objective_m.cc
new file mode 100644
index 0000000000000000000000000000000000000000..b3445753a375db9f59980a26dafee69dac276c43
--- /dev/null
+++ b/mex/sources/k_order_welfare/objective_m.cc
@@ -0,0 +1,143 @@
+/*
+ * 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/>.
+ */
+
+#include <algorithm>
+#include <cassert>
+
+#include "dynare_exception.hh"
+
+#include "objective_m.hh"
+
+ObjectiveMFile::ObjectiveMFile(const std::string &modName, int ntt_arg) :
+ ObjectiveAC(ntt_arg),
+ ObjectiveMFilename{modName + ".objective.static"}
+{
+}
+
+void
+ObjectiveMFile::unpackSparseMatrixAndCopyIntoTwoDMatData(mxArray *sparseMat, TwoDMatrix &tdm)
+{
+ int totalCols = mxGetN(sparseMat);
+ mwIndex *rowIdxVector = mxGetIr(sparseMat);
+ mwIndex *colIdxVector = mxGetJc(sparseMat);
+
+ assert(tdm.ncols() == 3);
+ /* Under MATLAB, the following check always holds at equality; under Octave,
+ there may be an inequality, because Octave diminishes nzmax if one gives
+ zeros in the values vector when calling sparse(). */
+ assert(tdm.nrows() >= mxGetNzmax(sparseMat));
+
+ double *ptr = mxGetPr(sparseMat);
+
+ int rind = 0;
+ int output_row = 0;
+
+ for (int i = 0; i < totalCols; i++)
+ for (int j = 0; j < static_cast<int>((colIdxVector[i+1]-colIdxVector[i])); j++, rind++)
+ {
+ tdm.get(output_row, 0) = rowIdxVector[rind] + 1;
+ tdm.get(output_row, 1) = i + 1;
+ tdm.get(output_row, 2) = ptr[rind];
+ output_row++;
+ }
+
+ /* If there are less elements than expected (that might happen if some
+ derivative is symbolically not zero but numerically zero at the evaluation
+ point), then fill in the matrix with empty entries, that will be
+ recognized as such by KordpDynare::populateDerivativesContainer() */
+ while (output_row < tdm.nrows())
+ {
+ tdm.get(output_row, 0) = 0;
+ tdm.get(output_row, 1) = 0;
+ tdm.get(output_row, 2) = 0;
+ output_row++;
+ }
+}
+
+void
+ObjectiveMFile::eval(const Vector &y, const Vector &x, const Vector &modParams,
+ Vector &residual, std::vector<TwoDMatrix> &md) noexcept(false)
+{
+ mxArray *T_m = mxCreateDoubleMatrix(ntt, 1, mxREAL);
+
+ mxArray *y_m = mxCreateDoubleMatrix(y.length(), 1, mxREAL);
+ std::copy_n(y.base(), y.length(), mxGetPr(y_m));
+
+ mxArray *x_m = mxCreateDoubleMatrix(1, x.length(), mxREAL);
+ std::copy_n(x.base(), x.length(), mxGetPr(x_m));
+
+ mxArray *params_m = mxCreateDoubleMatrix(modParams.length(), 1, mxREAL);
+ std::copy_n(modParams.base(), modParams.length(), mxGetPr(params_m));
+
+ mxArray *T_flag_m = mxCreateLogicalScalar(false);
+
+ {
+ // Compute temporary terms (for all orders)
+ std::string funcname = ObjectiveMFilename + "_g" + std::to_string(md.size()) + "_tt";
+ mxArray *plhs[1], *prhs[] = { T_m, y_m, x_m, params_m };
+
+ int retVal = mexCallMATLAB(1, plhs, 4, prhs, funcname.c_str());
+ if (retVal != 0)
+ throw DynareException(__FILE__, __LINE__, "Trouble calling " + funcname);
+
+ mxDestroyArray(T_m);
+ T_m = plhs[0];
+ }
+
+ {
+ // Compute residuals
+ std::string funcname = ObjectiveMFilename + "_resid";
+ mxArray *plhs[1], *prhs[] = { T_m, y_m, x_m, params_m, T_flag_m };
+
+ int retVal = mexCallMATLAB(1, plhs, 5, prhs, funcname.c_str());
+ if (retVal != 0)
+ throw DynareException(__FILE__, __LINE__, "Trouble calling " + funcname);
+
+ residual = Vector{plhs[0]};
+ mxDestroyArray(plhs[0]);
+ }
+
+ for (size_t i = 1; i <= md.size(); i++)
+ {
+ // Compute model derivatives
+ std::string funcname = ObjectiveMFilename + "_g" + std::to_string(i);
+ mxArray *plhs[1], *prhs[] = { T_m, y_m, x_m, params_m, T_flag_m };
+
+ int retVal = mexCallMATLAB(1, plhs, 5, prhs, funcname.c_str());
+ if (retVal != 0)
+ throw DynareException(__FILE__, __LINE__, "Trouble calling " + funcname);
+
+ if (i == 1)
+ {
+ assert(static_cast<int>(mxGetM(plhs[0])) == md[i-1].nrows());
+ assert(static_cast<int>(mxGetN(plhs[0])) == md[i-1].ncols());
+ std::copy_n(mxGetPr(plhs[0]), mxGetM(plhs[0])*mxGetN(plhs[0]), md[i-1].base());
+ }
+ else
+ unpackSparseMatrixAndCopyIntoTwoDMatData(plhs[0], md[i-1]);
+
+ mxDestroyArray(plhs[0]);
+ }
+
+ mxDestroyArray(T_m);
+ mxDestroyArray(y_m);
+ mxDestroyArray(x_m);
+ mxDestroyArray(params_m);
+ mxDestroyArray(T_flag_m);
+}
diff --git a/mex/sources/k_order_welfare/objective_m.hh b/mex/sources/k_order_welfare/objective_m.hh
new file mode 100644
index 0000000000000000000000000000000000000000..16eefb2db6bb94ef463001656b1ac6d8ce9cc470
--- /dev/null
+++ b/mex/sources/k_order_welfare/objective_m.hh
@@ -0,0 +1,43 @@
+/*
+ * 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/>.
+ */
+
+#ifndef _OBJECTIVE_M_HH
+#define _OBJECTIVE_M_HH
+
+#include "objective_abstract_class.hh"
+
+#include "mex.h"
+#include <dynmex.h>
+
+/**
+ * handles calls to <model>/+objective/static.m
+ *
+ **/
+class ObjectiveMFile : public ObjectiveAC
+{
+private:
+ const std::string ObjectiveMFilename;
+ static void unpackSparseMatrixAndCopyIntoTwoDMatData(mxArray *sparseMat, TwoDMatrix &tdm);
+public:
+ explicit ObjectiveMFile(const std::string &modName, int ntt_arg);
+ virtual ~ObjectiveMFile() = default;
+ void eval(const Vector &y, const Vector &x, const Vector ¶ms,
+ Vector &residual, std::vector<TwoDMatrix> &md) override;
+};
+#endif
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 38ab3118e066ade3dc2dcb08d9ced09c3672953a..395c9aa750c6b3d1881a2a898bc54e3a467b38ce 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -116,6 +116,8 @@ MODFILES = \
optimal_policy/mult_elimination_test.mod \
optimal_policy/neo_growth.mod \
optimal_policy/neo_growth_ramsey.mod \
+ optimal_policy/neo_growth_k_order.mod \
+ optimal_policy/neo_growth_ramsey_k_order.mod \
optimal_policy/Ramsey/ramsey_ex_initval.mod \
optimal_policy/Ramsey/ramsey_ex.mod \
optimal_policy/Ramsey/ramsey_ex_initval_AR2.mod \
@@ -597,6 +599,9 @@ optimal_policy/neo_growth_ramsey.o.trs: optimal_policy/neo_growth.o.trs
optimal_policy/neo_growth_ramsey_foresight.m.trs: optimal_policy/neo_growth_foresight.m.trs
optimal_policy/neo_growth_ramsey_foresight.o.trs: optimal_policy/neo_growth_foresight.o.trs
+optimal_policy/neo_growth_ramsey_k_order.m.trs: optimal_policy/neo_growth_k_order.m.trs
+optimal_policy/neo_growth_ramsey_k_order.o.trs: optimal_policy/neo_growth_k_order.o.trs
+
example1_use_dll.m.trs: example1.m.trs
example1_use_dll.o.trs: example1.o.trs
diff --git a/tests/optimal_policy/neo_growth_k_order.mod b/tests/optimal_policy/neo_growth_k_order.mod
new file mode 100644
index 0000000000000000000000000000000000000000..b61c9cb1a854637df7c4bd32c8b8d88a183913ef
--- /dev/null
+++ b/tests/optimal_policy/neo_growth_k_order.mod
@@ -0,0 +1,34 @@
+/*
+ * Copyright (C) 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/>.
+ */
+/*
+ * This file computes a second-order approximation of the neo-classical growth model.
+ * It is called by neo_growth_ramsey.mod to compare by-hand calculations of unconditional
+ * and conditional welfares and the output of the evaluate_planner_objective function.
+ */
+@#include "neo_growth_common.inc"
+
+shocks;
+var e;
+stderr 1;
+end;
+
+steady;
+resid;
+
+stoch_simul(order=6, irf=0);
diff --git a/tests/optimal_policy/neo_growth_ramsey_k_order.mod b/tests/optimal_policy/neo_growth_ramsey_k_order.mod
new file mode 100644
index 0000000000000000000000000000000000000000..e34db561368a93a3ad9ea86fedf3aed1999336f7
--- /dev/null
+++ b/tests/optimal_policy/neo_growth_ramsey_k_order.mod
@@ -0,0 +1,60 @@
+/*
+ * Copyright (C) 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/>.
+ */
+/*
+ * This file computes a kth-order approximation of the neo-classical growth model.
+ * It assesses the conditional welfare and the derivatives of the felicity and welfare functions computed by the k_order_welfare function
+ * and compares them to a by-hand assessment stemming from the results the model neo_growth_k_order.mod incur.
+ */
+
+@#include "neo_growth_ramsey_common.inc"
+
+shocks;
+var e;
+stderr 1;
+end;
+
+stoch_simul(order=6, irf=0);
+
+[condWelfare, U_dynpp, W_dynpp, U_dyn, W_dyn] = k_order_welfare(oo_.dr, M_, options_);
+
+if ~exist('neo_growth_k_order_results.mat','file');
+ error('neo_growth_k_order must be run first');
+end;
+
+oo = load('neo_growth_k_order_results','oo_');
+M = load('neo_growth_k_order_results','M_');
+options = load('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');
+end;