From e880d1bcc3b1153aebb83000e078f8a4d97746f2 Mon Sep 17 00:00:00 2001
From: Normann Rion <normann@dynare.org>
Date: Mon, 15 Feb 2021 17:21:43 +0100
Subject: [PATCH] Ref. #1680: 2nd-order welfare

---
 matlab/evaluate_planner_objective.m           | 225 +++++++++++-------
 matlab/ramsey_policy.m                        |   3 +-
 tests/Makefile.am                             |   6 +
 .../Ramsey/oo_ramsey_policy_initval.mat       | Bin 2896 -> 3654 bytes
 tests/optimal_policy/Ramsey/ramsey_ex.mod     |  12 +-
 .../Ramsey/ramsey_ex_initval.mod              |   2 +-
 tests/optimal_policy/neo_growth.mod           |  60 +++++
 tests/optimal_policy/neo_growth_ramsey.mod    |  85 +++++++
 8 files changed, 295 insertions(+), 98 deletions(-)
 create mode 100644 tests/optimal_policy/neo_growth.mod
 create mode 100644 tests/optimal_policy/neo_growth_ramsey.mod

diff --git a/matlab/evaluate_planner_objective.m b/matlab/evaluate_planner_objective.m
index a6df353260..26e23b4a15 100644
--- a/matlab/evaluate_planner_objective.m
+++ b/matlab/evaluate_planner_objective.m
@@ -1,17 +1,56 @@
-function planner_objective_value = evaluate_planner_objective(M,options,oo)
+function planner_objective_value = evaluate_planner_objective(M_,options_,oo_)
+
+% OUTPUT
+% Returns a vector containing first order or second-order approximations of 
+% - the unconditional expectation of the planner's objective function
+% - the conditional expectation of the planner's objective function starting from the non-stochastic steady state and allowing for future shocks
+% depending on the value of options_.order.
+
+% ALGORITHM
+% Welfare verifies
+% W(y_{t-1}, u_t, sigma) = U(h(y_{t-1}, u_t, sigma)) + beta E_t W(g(y_{t-1}, u_t, sigma), u_t, sigma)
+% where
+% - W is the welfare function
+% - U is the utility function
+% - y_{t-1} is the vector of state variables
+% - u_t is the vector of exogenous shocks scaled with sigma i.e. u_t = sigma e_t where e_t is the vector of exogenous shocks
+% - sigma is the perturbation parameter
+% - h is the policy function, providing controls x_t in function of information at time t i.e. (y_{t-1}, u_t, sigma)
+% - g is the transition function, providing next-period state variables in function of information at time t i.e. (y_{t-1}, u_t, sigma)
+% - beta is the planner's discount factor
+% - E_t is the expectation operator given information at time t i.e. (y_{t-1}, u_t, sigma)
+
+% The unconditional expectation of the planner's objective function verifies
+% E(W) = E(U)/(1-beta)
+% The conditional expectation of the planner's objective function given (y_{t-1}, u_t, sigma) coincides with the welfare function delineated above.
+
+% A first-order approximation of the utility function around the non-stochastic steady state (y_{t-1}, u_t, sigma) = (y, 0, 0) is
+% U(h(y_{t-1}, u_t, sigma)) = Ubar + U_x ( h_y yhat_{t-1} + h_u u_t )
+% Taking the unconditional expectation yields E(U) = Ubar and E(W) = Ubar/(1-beta)
+% As for conditional welfare, a first-order approximation leads to
+% W = Wbar + W_y yhat_{t-1} + W_u u_t
+% The approximated conditional expectation of the planner's objective function taking at the non-stochastic steady-state and allowing for future shocks thus verifies
+% W (y, 0, 1) = Wbar
+
+% Similarly, taking the unconditional expectation of a second-order approximation of utility around the non-stochastic steady state yields a second-order approximation of unconditional welfare
+% E(W) = (1 - beta)^{-1} ( Ubar + U_x h_y E(yhat) + 0.5 ( (U_xx h_y^2 + U_x h_yy) E(yhat^2) + (U_xx h_u^2 + U_x h_uu) E(u^2) + U_x h_ss )
+% where E(yhat), E(yhat^2) and E(u^2) can be derived from oo_.mean and oo_.var
+
+% As for conditional welfare, the second-order approximation of welfare around the non-stochastic steady state leads to
+% W(y_{t-1}, u_t, sigma) = Wbar + W_y yhat_{t-1} + W_u u_t + W_yu yhat_{t-1} ⊗ u_t + 0.5 ( W_yy yhat_{t-1}^2 + W_uu u_t^2 + W_ss )
+% The derivatives of W taken at the non-stochastic steady state can be computed as in Kamenik and Juillard (2004) "Solving Stochastic Dynamic Equilibrium Models: A k-Order Perturbation Approach".
+% The approximated conditional expectation of the planner's objective function starting from the non-stochastic steady-state and allowing for future shocks thus verifies
+% W(y,0,1) = Wbar + 0.5*Wss
 
-%function oo1 = evaluate_planner_objective(dr,M,oo,options)
-%  computes value of planner objective function
-%
 % INPUTS
-%   M:        (structure) model description
-%   options:  (structure) options
-%   oo:       (structure) output results
+%   M_:        (structure) model description
+%   options_:  (structure) options
+%   oo_:       (structure) output results
 %
 % SPECIAL REQUIREMENTS
 %   none
 
-% Copyright (C) 2007-2020 Dynare Team
+% Copyright (C) 2007-2021 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -28,88 +67,96 @@ function planner_objective_value = evaluate_planner_objective(M,options,oo)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-if options.order>1
-    fprintf('\nevaluate_planner_objective: order>1 not yet supported\n')
-    planner_objective_value = NaN;
-    return
-end
-dr = oo.dr;
-exo_nbr = M.exo_nbr;
-nstatic = M.nstatic;
-nspred = M.nspred;
-beta = get_optimal_policy_discount_factor(M.params, M.param_names);
-
-Gy = dr.ghx(nstatic+(1:nspred),:);
-Gu = dr.ghu(nstatic+(1:nspred),:);
-gy(dr.order_var,:) = dr.ghx;
-gu(dr.order_var,:) = dr.ghu;
-
-
-ys = oo.dr.ys;
-
-[U,Uy,Uyy] = feval([M.fname '.objective.static'],ys,zeros(1,exo_nbr), ...
-                   M.params);
-%second order terms
-Uyy = full(Uyy);
-
-Uyygygy = A_times_B_kronecker_C(Uyy,gy,gy);
-Uyygugu = A_times_B_kronecker_C(Uyy,gu,gu);
-Uyygygu = A_times_B_kronecker_C(Uyy,gy,gu);
-
-Wbar =U/(1-beta); %steady state welfare
-Wy = Uy*gy/(eye(nspred)-beta*Gy);
-Wu = Uy*gu+beta*Wy*Gu;
-% Wyy = Uyygygy/(eye(nspred*nspred)-beta*kron(Gy,Gy)); %solve Wyy=Uyy*kron(gy,gy)+beta*Wyy*kron(Gy,Gy)
-if isempty(options.qz_criterium)
-    options.qz_criterium = 1+1e-6;
-end
-%solve Lyapunuv equation Wyy=gy'*Uyy*gy+beta*Gy'Wyy*Gy
-Wyy = reshape(lyapunov_symm(sqrt(beta)*Gy',reshape(Uyygygy,nspred,nspred),options.lyapunov_fixed_point_tol,options.qz_criterium,options.lyapunov_complex_threshold, 3, options.debug),1,nspred*nspred);
-
-Wyygugu = A_times_B_kronecker_C(Wyy,Gu,Gu);
-Wyygygu = A_times_B_kronecker_C(Wyy,Gy,Gu);
-Wuu = Uyygugu+beta*Wyygugu;
-Wyu = Uyygygu+beta*Wyygygu;
-Wss = beta*Wuu*M.Sigma_e(:)/(1-beta); % at period 0, we are in steady state, so the deviation term only starts in period 1, thus the beta in front
-
-% initialize yhat1 at the steady state
-yhat1 = oo.steady_state;
-if options.ramsey_policy
-    % initialize le Lagrange multipliers to 0 in yhat2
-    yhat2 = zeros(M.endo_nbr,1);
-    yhat2(1:M.orig_endo_nbr) = oo.steady_state(1:M.orig_endo_nbr);
-end
-if ~isempty(M.endo_histval)
-    % initialize endogenous state variable to histval if necessary
-    yhat1(1:M.orig_endo_nbr) = M.endo_histval(1:M.orig_endo_nbr);
-    if options.ramsey_policy
-        yhat2(1:M.orig_endo_nbr) = M.endo_histval(1:M.orig_endo_nbr);
+dr = oo_.dr;
+
+exo_nbr = M_.exo_nbr;
+nstatic = M_.nstatic;
+nspred = M_.nspred;
+beta = get_optimal_policy_discount_factor(M_.params, M_.param_names);
+
+ys = oo_.dr.ys;
+planner_objective_value = zeros(2,1);
+if options_.order == 1
+    [U] = feval([M_.fname '.objective.static'],ys,zeros(1,exo_nbr), M_.params);
+    planner_objective_value(1) = U/(1-beta);
+    planner_objective_value(2) = U/(1-beta);  
+elseif options_.order == 2
+    [U,Uy,Uyy] = feval([M_.fname '.objective.static'],ys,zeros(1,exo_nbr), M_.params);
+
+    Gy = dr.ghx(nstatic+(1:nspred),:);
+    Gu = dr.ghu(nstatic+(1:nspred),:);
+    Gyy = dr.ghxx(nstatic+(1:nspred),:);
+    Gyu = dr.ghxu(nstatic+(1:nspred),:);
+    Guu = dr.ghuu(nstatic+(1:nspred),:);
+    Gss = dr.ghs2(nstatic+(1:nspred),:);
+
+    gy(dr.order_var,:) = dr.ghx;
+    gu(dr.order_var,:) = dr.ghu;
+    gyy(dr.order_var,:) = dr.ghxx;
+    gyu(dr.order_var,:) = dr.ghxu;
+    guu(dr.order_var,:) = dr.ghuu;
+    gss(dr.order_var,:) = dr.ghs2;
+
+    Uyy = full(Uyy);
+
+    Uyygygy = A_times_B_kronecker_C(Uyy,gy,gy);
+    Uyygugu = A_times_B_kronecker_C(Uyy,gu,gu);
+
+    %% Unconditional welfare
+
+    old_noprint = options_.noprint;
+
+    if ~old_noprint
+        options_.noprint = 1;
     end
-end
-yhat1 = yhat1(dr.order_var(nstatic+(1:nspred)),1)-dr.ys(dr.order_var(nstatic+(1:nspred)));
-u = oo.exo_simul(1,:)';
-
-Wyyyhatyhat1 = A_times_B_kronecker_C(Wyy,yhat1,yhat1);
-Wuuuu = A_times_B_kronecker_C(Wuu,u,u);
-Wyuyhatu1 = A_times_B_kronecker_C(Wyu,yhat1,u);
-planner_objective_value(1) = Wbar+Wy*yhat1+Wu*u+Wyuyhatu1 ...
-    + 0.5*(Wyyyhatyhat1 + Wuuuu+Wss);
-if options.ramsey_policy
-    yhat2 = yhat2(dr.order_var(nstatic+(1:nspred)),1)-dr.ys(dr.order_var(nstatic+(1:nspred)));
-    Wyyyhatyhat2 = A_times_B_kronecker_C(Wyy,yhat2,yhat2);
-    Wyuyhatu2 = A_times_B_kronecker_C(Wyu,yhat2,u);
-    planner_objective_value(2) = Wbar+Wy*yhat2+Wu*u+Wyuyhatu2 ...
-        + 0.5*(Wyyyhatyhat2 + Wuuuu+Wss);
-end
+    var_list = M_.endo_names(dr.order_var(nstatic+(1:nspred)));
+    [info, oo_, options_] = stoch_simul(M_, options_, oo_, var_list); %get decision rules and moments
+    if ~old_noprint
+        options_.noprint = 0;
+    end
+
+    oo_.mean(isnan(oo_.mean)) = options_.huge_number;
+    oo_.var(isnan(oo_.var)) = options_.huge_number;
+
+    Ey = oo_.mean;
+    Eyhat = Ey - ys(dr.order_var(nstatic+(1:nspred)));
 
-if ~options.noprint
-    fprintf('\nApproximated value of planner objective function\n')
-    if options.ramsey_policy
-        fprintf('    - with initial Lagrange multipliers set to 0: %10.8f\n', ...
-                planner_objective_value(2))
-        fprintf('    - with initial Lagrange multipliers set to steady state: %10.8f\n\n', ...
-                planner_objective_value(1))
-    elseif options.discretionary_policy
-        fprintf('with discretionary policy: %10.8f\n\n',planner_objective_value(1))
+    var_corr = Eyhat*Eyhat';
+    Eyhatyhat = oo_.var(:) + var_corr(:);
+    Euu = M_.Sigma_e(:);
+
+    EU = U + Uy*gy*Eyhat + 0.5*((Uyygygy + Uy*gyy)*Eyhatyhat + (Uyygugu + Uy*guu)*Euu + Uy*gss);
+    EW = EU/(1-beta);
+
+    %% Conditional welfare starting from the non-stochastic steady-state
+
+    Wbar = U/(1-beta);
+    Wy = Uy*gy/(eye(nspred)-beta*Gy);
+
+    if isempty(options_.qz_criterium)
+        options_.qz_criterium = 1+1e-6;
+    end
+    %solve Lyapunuv equation Wyy=gy'*Uyy*gy+Uy*gyy+beta*Wy*Gyy+beta*Gy'Wyy*Gy
+    Wyy = reshape(lyapunov_symm(sqrt(beta)*Gy',reshape(Uyygygy + Uy*gyy + beta*Wy*Gyy,nspred,nspred),options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold, 3, options_.debug),1,nspred*nspred);
+    Wyygugu = A_times_B_kronecker_C(Wyy,Gu,Gu);
+    Wuu = Uyygugu + Uy*guu + beta*(Wyygugu + Wy*Guu);
+    Wss = (Uy*gss + beta*(Wy*Gss + Wuu*M_.Sigma_e(:)))/(1-beta);
+    W = Wbar + 0.5*Wss;
+
+    planner_objective_value(1) = EW;
+    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;
+    planner_objective_value(2) = NaN;
+    return
+end
+if ~options_.noprint
+    if options_.ramsey_policy
+        fprintf('\nApproximated value of unconditional welfare:  %10.8f\n', planner_objective_value(1))
+        fprintf('\nApproximated value of conditional welfare:  %10.8f\n', planner_objective_value(2))
+    elseif options_.discretionary_policy
+        fprintf('\nApproximated value of unconditional welfare with discretionary policy: %10.8f\n\n', EW)
     end
 end
diff --git a/matlab/ramsey_policy.m b/matlab/ramsey_policy.m
index 720e2d8337..69b2a03b60 100644
--- a/matlab/ramsey_policy.m
+++ b/matlab/ramsey_policy.m
@@ -1,6 +1,6 @@
 function info = ramsey_policy(var_list)
 
-% Copyright (C) 2007-2019 Dynare Team
+% Copyright (C) 2007-2021 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -21,7 +21,6 @@ global options_ oo_ M_
 
 options_.ramsey_policy = 1;
 oldoptions = options_;
-options_.order = 1;
 
 %test whether specification matches
 inst_nbr = size(options_.instruments,1);
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 46f29f5367..07dcebb9b6 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -109,6 +109,8 @@ MODFILES = \
 	optimal_policy/nk_ramsey_expectation.mod \
 	optimal_policy/nk_ramsey_expectation_a.mod \
 	optimal_policy/mult_elimination_test.mod \
+	optimal_policy/neo_growth.mod \
+	optimal_policy/neo_growth_ramsey.mod \
 	optimal_policy/Ramsey/ramsey_ex_initval.mod \
 	optimal_policy/Ramsey/ramsey_ex.mod \
 	optimal_policy/Ramsey/ramsey_ex_initval_AR2.mod \
@@ -573,6 +575,10 @@ XFAIL_MODFILES = ramst_xfail.mod \
 MFILES = histval_initval_file/ramst_initval_file_data.m
 
 # Dependencies
+
+optimal_policy/neo_growth_ramsey.m.trs: optimal_policy/neo_growth.m.trs
+optimal_policy/neo_growth_ramsey.o.trs: optimal_policy/neo_growth.o.trs
+
 example1_use_dll.m.trs: example1.m.trs
 example1_use_dll.o.trs: example1.o.trs
 
diff --git a/tests/optimal_policy/Ramsey/oo_ramsey_policy_initval.mat b/tests/optimal_policy/Ramsey/oo_ramsey_policy_initval.mat
index ff829dd9dc5fbf2a4d9b0d649c50e1537c7a672f..6c363f3b29c47f9f1450d7e640b175fd97661162 100644
GIT binary patch
literal 3654
zcmeZu4DoSvQZUssQ1EpO(M`+DN!3vZ$Vn_o%P-2cQgHY2i*PhE(NS<NN=+<DO;O0t
zvr-68O;K=5O;Rv4S1>fPGBvX@HC8Y(Ffvpi5}e~<fB;5MUw#G#1~(oChKf0FA~Qoo
zj*G2}lu+sFT69p=W${Bs{{Z1nCQ88)(^a&(ypJf{GBG}T!GQZ4k8-j@yToK0-x|4o
ziSin`#XFdeJ>WjHWr5J<6&!jy9t0OU1|E^HZE%Y(_w_e6VeWY%^Z9Y*`ISGP-COnh
z?yl{tRxVm@%&;Jffx(E)+Q64*&7(svm3BWnYBHB^aY64qgTziN)0GB^OQxAf6mI37
zu&Ayq#_aP#lgz@d13fHD)k~uea8F!w_{$PCPu^Kv3q!5m1Oy#es36DOsO`)1fbGz|
zbt*h(e^u;OU@@=t@efd0n7?K9B#-t&*`p~(JA!A-x+2~3ICO_wbMmq@8!4viX2wyK
zGRIQisH_VW@Dbvgs}{Iiop;5$ib;;kF9e>Q<DnJyzC%5*H+w5f>&v7^oUM%}M&2!}
zCv#kw`I4Kp_qW5T-tG6<=USvJ)n3ld%E#?-$>?cE;Nq$cHqLJ<Qkdj8Gw%4W>6o}Q
z*M7FiFUEWMThHwLo4wV0BLC`iBi)53y(UE|FIzmtr}jo^#-C$dD?M#P3s3rPQL*L?
z{8Y!XGRf=AcTG><#~zD1mu*_k^4|QW$(uQc9Z!EPWW8p(dA)!9wR;DBl}aaist6e>
z>luZwi&vTEk-KN%+7G)JA67-?u9I~US~<J>Or7kl!^(w|cNxSa&6M^3RQzGXFLod9
zcZDlH-P~MlkbZjMnsB)*4xirtk2;q0YY(I3jCkEGR}0fimwR^i^IKZHihKE{;{@-t
z+qU2B{yBb@t%<Pu_9-q>*QHc!t6aQR$Nm1mW3AyqciG*e-L4#8!xQnl*R{6z&b#DU
z6+3nFw(pCpo|3qzsb#a^(`9p1r*W=~TP)2rW!j^79%~!P{d@juyUgl0KV`4O$Gc;<
zV%WTE*(Z9he>E-HvSrJMC!c&~b)+tBJ9$`U%H&<wm>f=iko*55$S&{o(K8R<B~I)8
zzVB|lWVN`SZ{~|h9qG}#bpF&_{BlXgBFXPZ_Ur#USnJj5IZVFkEy(?Bxo_ROk8a;H
zKh>5W`7Ecco$%pF`t8H#pB>u6+I+75$dA-NyR6oJ`5kr3=8SU2y4z~I9>4zly|k{F
z%l3k8dECh-TWf5C-KDS2+xz$FvE6p;kMG@J_uf@w`~T9k*e!{dA8Ly{zV1*H{pD|k
z)4hA8cV88Lb9q)cb<gSZc4}$0FC5-4TB)v=vwZF2tcBTLepSb=mWpvJa?hOc$-mO=
z$%7-g+Uk*Ar_=INuixKLGw1v4Pcq)o_3<taKmK?y?lzu%cDHk|eK`APxpQ0GZBH4k
z-mW{L{TpxfilkSYtag+??OGn+)ARJ&$Jff?%>6N4a#LqZ$nLQ_tv}^g^z)g2RlhHs
z@pwz8@#nl>F`WT&Z$C~hs5>3~adN}6-}bXIb>BXp_ATA}_u8$0C*M7_X9?f4)U>GR
z^NG>t!=umpMxS@CzIJeC>GsaWpEvBSOS87;KWI@F-r76QruXX2zbl{U`K-CmURL+6
zslv#yZq2*Yzfp=}PxK!>6V$8Q_j!-~a&hC{Lt7;lnop3rbK%Lolib!0=kLsXBopVn
zlaVXN_gF$zuU~qHa!a{(nS+w!QKRqbI$`dDRY9WZ0(=*azZBs-bb8ibhOPrw&zN+I
zxp7^;=zFJ1HtLr7xgVXkUS7XfQ9t#5eC>LUt-HfJSU$$Gd~B_>+STv<=2&^^7MD3C
zr~I!R@n;e7J0|?g^49WmKLX!=Ic&M#BJy;p$h_lkepp@jv}<?6v0IxDt@|B#L8H-K
zzBMPqA|m1PU#92ZX0<GBN%_hD=S+QM;=zP<|EKalzJKn!9JBva_r@42<0WhIzV6y0
z+`u>go~`cR8xnQ%rxg5kQ23p_uH<Lt)o%jveg9RS=fA!`#~{S;+tr$%2SjIm4|}Y1
zOYHu!`KN#0{<T+Mq;LM`mk~dne_9`V;-N94*=BQ=*wCo+{@MSxNprUuTFD-sznbrF
ztIeL;)3-bNw|%XSaj)4RZM8!B?V%mkJ9l{Lz3V<zANuNV?ECnghtDf(R~~xa;}<{U
z{m$d(-j)0_4ib=`pMJ>Z2>bQf0`mFs${o9>cSj#QS{wLt_p|EfITu1df8HSy|3h@~
z<;m9UZ?<l>p8Nl&&5`5JXSyoc-8e7%@#8nYrXACF?g-woLtpEEx$W^IFL&L0_I%-N
zm;Tg}s*Wt4%=B`r|8aK(<>x<NVsM1rHc3!EpI_|A+|x&EySslc`n3OE`Rx;v7HQSV
zP5!L^EwSX+-7i}e%gUC2w!g~$@ndD1oI}lpa@9wV--b5rxV&S>!I&7c{||GY=-22s
z=B#0O{nPNX<)53selk6O&#Qla`M0(GlSJ3m?bFk}&+gLnYp)G|{JeJk^M|KDfB64%
z<&3w#YTk$Eyf4@@^VpyMPE|(MWB=|gbe%8q>xkiNeFdHRB<=U-Bjw^9Yn?4qJ%7~S
zEIxCmkZHq)7ex!@UT*&LN}%h#jZ;MZ6Gf}{ZYh8JF6{Wqc<7ew{)PYS`R;5_s?<Hi
zA?Wl$VIAAUeW$EBRzD9@KJ%hX`r+2Yp&Y08@8~Rue!|<*?H=d8qx^ztw#|R97Up|9
zyDk(z5^dSNb-G|0zk6ka-)(JyV{9=lGs=J1IO)H&Y*}4%TF}kbBKd0nT(=*Y?=)Ks
zcb^uhn-<&^5LY6Ucy_5!i+@VknJe=88Rd85dF?Gt(iK@ax1X!Ona*gke_7b?b-s<y
z5@+|7N`2vExh(H`P)L}iw7br7uZ85S`DZ<L+b7DgEZiOQ@}C=rQ_yqulamxUK8T-P
zaozbrpy;;E^P18g|6;U1sll;r&aO)v@5gW)@vfTI@7dxI^-R9Ro7?fI-2J7SO$B^x
zx2GNd%EIAxR%&)$U;6>=(%2LR&Vx7ayt7>$Dxh;FI`H%S2o5%0-S->28WjG#XqLY`
zrKdLI;SYZifhn>b`)7+WH-1}tKjv~FM_%eh-utf38#Fy+r}jI3czxNFx%Qbu!iPg{
zd7VuO#Rm`el&3igRDKRPva???A~ZB7%1pTF{$j1#W$q3Wew!q%ywAhYHUC=pj(V1Z
zGg)^V8qA#;{&K5_0dJE!qaE{|<xUeO{pnqPUgU@Rx7|C{cJbZayhm-CqK?)Yhl<1Z
zyQ1s&G#Ip6^c$@IZP_yI(q<0tf3<82xBoObbo*5>$K=p3JN;PALyg^jN;R8yvLs&r
zs=nrNOn>_yeg=mB|C19E5+-EWCvz{6u=AM5c3I}Pf78oZmxL3n4nMP)@lUyg$4C5+
z1=F&~dopbo^7$KID(~lOddaFJFzeteh8gn;Q{`v;d$?2FhdbsT@3QzOZ<#I&r<65i
z_TIjI@Jr{^+Xue1eOsFlwJ?Qy2H%O>JcjE_ChN+dknQ2&$x--j<d7kd(E3Tz@tI?$
zzyp>fc7fvxw)M4BuhzBke*g7tW%RrB9S8oreh~cmFAr;~^t)U$=B{v=lGz!%<~ug3
z-d=N!MI&zC%i@`a6Yl^0l(Y6?U0X|0@y*V*-;!cBf#tUn7)1WZKi?W!Z-2N^`qs`*
zRlg;rrtRNj8qU5V@6z|W+P8ms_?P+5sfsMNZ+31FPhQRz@UQ-yxP7U@e)ZaIT~ZVN
zAC8(HEV|;|@&54rQyZVLM4tK|$aHFbgyw0^CGto73xk60oIfBC`{$%DqsF`g>3_ob
zR-F4-8L>zwAm-78*WQ^6EZNri=mt4_jAc)k_F#IcDqMQMllkSpTh+UdHk53=)-1Wm
z;cxq|Pf;!nuOdI6yUxY@(>jlR8wcar%K97SyA&p@k6F6Mm*GVDyx6h^`4xGSrKgL@
zPq26Ve7I=d-zRO=Yo2}g-sj2tO<gkgRH^9ew|mpL|K7G+d2L_lkK>xQol@0PO>b^j
z%U!W4e9hFFpN6Y0)tP?2#Jw^}^xM_;=X(U}LelQvSayBg=IvMh-KgijAAZ{F&GK(Y
z=FG8Je`DYF(#MyZj=Fx!-f>?1{_UWQ_mSm=)1NM7s(yO(RZ{rAXSd~}_m$k8H^Zjd
zc5b-0YW|ZO{w3$GRn3}~Yyb7;WjE74o$7fv(}nHdo!U_U$xTgRtKQu@&#5NoZyhlC
zeN1=u|G%bwTgCJK?kn8o)%ABz|I>f%QzyOs?&7f`?waqW?|cu90*ZXM#LwPYv9`di
z>2G2Gvy{-e6~FtgzmJ~s{MfD4Td&rM1a@tFzNhlw@9XC@e&64$Q*ZClx8;0p@w*BY
z@x=Gik!MzX(mG=o|GwhaRpo`>e5Wtdlw8NX(LW|UDs->Ojep;IY=eAOH@%tp=u(>3
zd%i8~t-G0--xZ`cwOJPJ;ASX1vv}=R<Fh*a2QnD)Z@=9G>P-|df_f86L-zU|HW0Ae
z?KwqBs<U@uW8=h)%0fa-3zHHSM2oIuQ~D@8iO0rCzVrA4)(H}ll`QAJ?c6s@CG_x0
z@$cvEmzUMV&tW`Z!gxS_$&2I-V!=9JX50vzy+q>3RiERMzntIv+|>I1&NF3Wbq-_e
z{rd__`y%EtAL{-pA-Sec{e}3^#k<yi+7xm%$??zn6QBQ9nU!Vl+H=g`^?UX5Bw>9`
zbEP*jwXXjf|2cm?H*M=GuF%6q`9&4y=8KrT@^zn4t(l$isk%mF#{Sh;w71OSc(g<I
z+G)Wp<&sZQxvxFE+-I35DeZr7e)6Gjvo_6IcW=hWa}xRX4|)}@^d4{d{bSnQTK?BF
z_Brp5V>J6J_38f6yefYUJ2&&H4NV4J49hq-tcqB_`EATkyPrGfIef2d*`xTwyKO#8
z#rJi5>mNJ+{vN!^a`{Qc*ZHsGEB@90+OB;6u59&-XW6Fr*RD6;;_qCk{bK+6=RNA0
T&%T?!zhD0@J9oe2KMNB8LbxXo

literal 2896
zcmeZu4DoSvQZUssQ1EpO(M`+DN!3vZ$Vn_o%P-2cQV4Jk_w+L}(NS<NN=+<DO;O0t
zvr=#?%2e<w%~3G4P%t#MGBmO>v`{cIFf>sh5-`93qo*%F0|UbeE(V5*Id5lhhD<&#
zc3i#1HDm>sR@a;E28WG{nw;9AE{3%(V$xcm=+NNcVW6m}kPvWiTh#XO0LIm!DIYtT
zwASkUjEIdAXUP?u`aCW&W#u;3EuYjT3mC51`n=+t^?A$mPuAbho%w%m?#%bnKdvw`
z%-zMHq1@*b%IMKG>5AIdl~F4;-3;N=?eXqk)V6sh&wmZA=A#`CF3<3C5ZvDV`{qLV
zT@AiZjSf}#d`MfknTaij<3R%>_qEOimzfp`WC(uQ_m82<Wxv_Xc)cSIQjdgOxVaYx
zF8+0tck$&2=Ckeo(iMwL;#h9aKl}3fi3co2?62g{E!h}R!Z-D;xYPIkgby!6wshE7
z1-R#$J@{Iu;r;vjJ9f!E_xK9~56Bt3*>UevmCS@*>E%Y7{!afP+BomZHwXJ1KZa9D
z8$ua*#O7PMOuta~`=eXug9{NcZY<4b7hS4gVw>#v*rl^o!o*%}F83m~n8-*@E-{|F
zXQWq2p3nJk<VBCdmc-_6en}~vcl8gRbXCO9Iq=~3O0Fm88;`h$g@y!4x&{UOYq#=f
zUbAA=k~`n$?BHG4+$Ufqng7b?u+^#*SJhObnQI!lCT;54GIxnX*qf`hw~MuZNLO~g
z^$(nDzVVTxzrVE2^ti&JFPp<8%#)w<8J<<{<qJA{J28LHoM#)?afau3=`TrL7yY!<
z;fa;u|D!iI1b^t4dcOPTskZ{_*)w1K=qmd7b>hagr<zxPT)19C{j$YLt<?1^;=JuY
z7e#O^HGlVZQt5K07!TF(+h6kBp4O@^`v0ZIvD~xa#k^Bba(qqn%5%0j-I$fHDO|&x
zW@9Vh(K#>ab434buNVERHFs~*?q0y&<*j|n*}z)g%$22~;NKC;_yE4&lFx6f+}W`&
z+xpkv<^#(EVwd{b&Y7nDsmJ@yo4CpBhp+E?x#aIvj?2yGqdqU6Cb0K&q48Tq`H*PO
zlNn(hHG&1~YP-Ir?EV|X@b<CJYugJJi@kq8%09P6cv|uC+v~z@U2?*^?AGROnp3ed
zP9XAcWmk@^DetBovx0Rs*6J=>da!MI;#Kx|ck?Ck>L&YE?>X-HxTfC2qBxCn?x$_P
zvtI9AVfETe!1b;BWsNQNE?ucVyH;=C)ERF3=n`9c{@n>RAL>?loSpP)ht4zauEU0k
z>Yt^<FTTCCRr1ka{&{Julbucq^OnnbI{cfX&Gb?Ddj7q|_hgnE9oX|WpfpkUM8>(E
z#UYuy|H!pmJNo|C^oeX{$7Pw8{(iCjZ)a&}^ZVo4*MFb$x6XX0GGpDoTMc&@MVHLl
z>nc-!UEs4Ce}P<WpzNC(-m*Bh+y57SGcNC_)Z{bJS6p@2R{ZeX?+?#?EO=h5^8CTs
z?`odg`160;{f@8YpQX9K&wo<)dL8ljJNt_}erkDU-~So+V)B2rf9W3_$}>IJzco}2
zzCBgy@r#XB!A(c|HcisKRJQiKQ<~^&y_Sh#mv^5Ri0q3^sL%Qn>QJ#e`-pmMoW|*X
z{ap5<*3CWZs|+n~UkqM#-+JHnzd?_D>YtxHIz?79k^9?4510HJ`IUCzR&~p_+<g*u
zb5iN-IS=hm@2_1v#kT0p^Y>f!$}3drsk}Yd|M;2GY>o-G`pkXPuTT8hf96n?OaAY#
z?GgN057<9!-S4fZ7p(Vg_bGeJZI6FEs@C09`GLh-c%rX~OZL1MkL3&hy)%sPpSHYI
z<y!N7w?x~S_L(zwyx01F@#wBf_8Vt6Ri?B}$#pl|7;{%^e}DavpV!5cY$NB@uw83z
z%Zl)SCa)}R?krW_YVvOPsrr4-4kdisp?1^eO$=vxity$pF{58U+~@3{#&&J;)0Iz~
zzCLJiSfl>=_Q|7Lrytuo^;m8APy2b(UIxVeb1~!J`{?$q6KXP(f@FUGaGz7Z)(oU4
ztmW&1pXL$%&!+dU))x<-^Z!-ir|J9GpE+>i{p#=Exuf)#9(<oXDdyk%_lv)Ozxez2
zi@&vB{@2Jq{de~L{JsA#zgQ&YmmE-eV7}JBX(tZtm-*wp*lV9&*Vn`CKm1g#8BUIl
zS6}&Y>h1a?-zOBD2~rRWei<gdZSPs@kAGeAelp*V4c(tyeedM+8Gm}xXS_XLcld3+
zXsJiU^}d~5Eu8aMw<pc3yU=uV&pye|4_egDv`_Unc)jgk*UXn+C5#U~6inOq)Y;Rf
zn|+4c=J&~2db`_a9^TbGGrRpWL-OYL32ae!kLj#gbNouu-6K9`s}Fyv{I>bXz8dD@
z;De<;`{zWR(XBL{G38VCjJ3R1XXsB=J@Y6|v+{N6{5NH}Y4_YXoTOjxztGq0B&&8=
zN0p81x!f+BYemc~@9KXXc%{K%W&XcJ*^!~ucHT!>Z)L>-t94Ei9~29EdY@^2R@ksf
zp|Y${QR4rRr}txJG@dM%j176_Frn`6J&C_x1vF|*SLP*Juykz|+t#bY@otX6?QLc}
zO`_)C<*P-w4j#XE&hwA}$Cjkvy=x={{sb*BTxrSC)Z1DeRnO6KpeXI`krXSIJn?yR
z`}rOI%}n(+`{5iAEb__Z9zWAH_x&qXE-GZ$c$Uon)cRn8?&ak(xDGB2|K00*-BDnn
zYpKWg#s>wSFAv_lCs6S>YU8VCQXG@!&*{FkQL5>=9ly`raL0fTSy#=9{W;F<eYI=x
zGsc7G3oSz)^gRpxe1UsS%Y&s;S5Cae#gW8PZR$JyaDBnO*0uZQ{&=slHg)Tar>W%z
z^3fhM3jc|v@$8k`9Pn6e*@|c4XI{OpJ9b6%-hrH?xZ@%$o?o`yFr9I1_WHzqHZi@H
zdT&=PnSV;ubN|}u2Dg|0X|e3uq0xDIcZ1}cf;gkPLzZVgob&QII*p5S{<P$<vgwI(
zm!=4u(LQ>oJhtWL&&W=dx?I6$!eNTf)@7y}e$B7#uT<R4U|jxMveW*N$ujBJQ(XF|
zzyERV^*fOxc9ORJ3`=>fT`s@nESP>L_FIM?<Ds8Bd6z#EPzVp$mhnA@>4kRNxz4Ss
z4F*R4CF{jDZ8weU``f{=Z(?z2*?E-%3)?S?Ke;7vV5NHA$(j5U{qFV{*k(UDpR>k3
zq4U?<pT*DHjVIK<f43y6_{#ZQ&X?ZCORgWTd1obi^LBGkukXaWmiN7`vV490CZ{T~
zWyhho_e&MC&Y!v#=(@-IZ(8F*#;*5g|K?7#`}BL&`Kndj&;Ea0KhH!qC3~v1d6u%Q
zbb709@I~8zkZD=QtXr4A-<CDEXNU2&g&9JR*t&QBR$KJ%sdYrONnpL|wW|@lDgO(y
ze@^4<$Q9MAKCUd~DlOXMY&Cz$CvSU)zn{c=w*|O=)BF7G%*G?P*X`c^bH1rwol}Hs
zQO6de8yw=x{VKlxofomgbMa^Q?6ndZ1vXpi{PRq%^v=C><Dn?)<`eF#UruHT2$*<n
zQQ>_4$3knjM@{>XklH`jJJ^lo&-JqT-MjZcJ(SB?*?*|N;K4@2@SVKNuOB{Ic%AEQ
z_MwIMCcd1!?*F;Ar!(FCMSgEs(3*PrqGqYjkFV+r#rm0?AI>}8|H#+9^}xlI(=%mT
z4}PdSZDL*(BVZ%5*48Y&QsGX;uPdwfr(Jybqs9L7H!1c5Tl!9%*XH?~_%KU%a+J3H
zjw7{OdKVh=%B9PkyvXsbHFSR7H}j244y`?#8{NaVrd6$TqwK6UyMilcUoE|xzBi$~
z>u>k&f7uGgA1CX-d%dph^tA%dW4`kyq}JNKa9hgWv)<EePURulU-KTg=3DF%y7pds
z-l~EeR+Y8fch<+MRGznJTzEa!G4Jm2g{vk{+gHeUbf3vZ+w12UKJJfIttkfpblR;U

diff --git a/tests/optimal_policy/Ramsey/ramsey_ex.mod b/tests/optimal_policy/Ramsey/ramsey_ex.mod
index d28882f970..04e1244626 100644
--- a/tests/optimal_policy/Ramsey/ramsey_ex.mod
+++ b/tests/optimal_policy/Ramsey/ramsey_ex.mod
@@ -53,11 +53,11 @@ end
 
 load oo_ramsey_policy_initval;
 
-if max(abs((oo_ramsey_policy_initval.steady_state-oo_.steady_state)))>1e-5 ...
-    || max(abs(oo_ramsey_policy_initval.dr.ys-oo_.dr.ys))>1e-5 ...
-    || max(max(abs(oo_ramsey_policy_initval.dr.ghx-oo_.dr.ghx)))>1e-5 ...
-    || max(max(abs(oo_ramsey_policy_initval.dr.ghu-oo_.dr.ghu)))>1e-5 ...
-    || max(max(abs(oo_ramsey_policy_initval.dr.Gy-oo_.dr.Gy)))>1e-5 ...
-    || max(abs((oo_ramsey_policy_initval.planner_objective_value-oo_.planner_objective_value)))>1e-5
+if any( [ max(abs((oo_ramsey_policy_initval.steady_state-oo_.steady_state)))>1e-5, ...
+    max(abs(oo_ramsey_policy_initval.dr.ys-oo_.dr.ys))>1e-5, ...
+    max(max(abs(oo_ramsey_policy_initval.dr.ghx-oo_.dr.ghx)))>1e-5, ...
+    max(max(abs(oo_ramsey_policy_initval.dr.ghu-oo_.dr.ghu)))>1e-5, ...
+    max(max(abs(oo_ramsey_policy_initval.dr.Gy-oo_.dr.Gy)))>1e-5, ...
+    max(abs((oo_ramsey_policy_initval.planner_objective_value-oo_.planner_objective_value)))>1e-5 ] )
     error('Initval and steady state file yield different results')
 end
diff --git a/tests/optimal_policy/Ramsey/ramsey_ex_initval.mod b/tests/optimal_policy/Ramsey/ramsey_ex_initval.mod
index 170b956f99..9f0d8e1326 100644
--- a/tests/optimal_policy/Ramsey/ramsey_ex_initval.mod
+++ b/tests/optimal_policy/Ramsey/ramsey_ex_initval.mod
@@ -40,7 +40,7 @@ end;
 planner_objective(ln(c)-phi*((n^(1+gamma))/(1+gamma)));
 ramsey_policy(planner_discount=0.99,order=1,instruments=(r));
 oo_ramsey_policy_initval=oo_;
-%save oo_ramsey_policy_initval oo_ramsey_policy_initval
+save oo_ramsey_policy_initval oo_ramsey_policy_initval;
 stoch_simul(periods=0, order=1, irf=25, nograph);
 if max(abs((oo_ramsey_policy_initval.steady_state-oo_.steady_state)))>1e-5 ...
     || max(abs(oo_ramsey_policy_initval.dr.ys-oo_.dr.ys))>1e-5 ...
diff --git a/tests/optimal_policy/neo_growth.mod b/tests/optimal_policy/neo_growth.mod
new file mode 100644
index 0000000000..2c5165c8c7
--- /dev/null
+++ b/tests/optimal_policy/neo_growth.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 <http://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.
+ */
+var U k z c W;
+
+varexo e;
+
+parameters beta gamma alpha delta rho s;
+
+beta = 0.987;
+gamma = 1;
+delta = 0.012;
+alpha = 0.4;
+rho = 0.95;
+s = 0.007;
+
+model;
+c^(-gamma)=beta*c(+1)^(-gamma)*(alpha*exp(z(+1))*k^(alpha-1)+1-delta);
+W = U + beta*W(+1);
+k=exp(z)*k(-1)^(alpha)-c+(1-delta)*k(-1);
+z=rho*z(-1)+s*e;
+U=ln(c);
+end;
+
+steady_state_model;
+k = ((1/beta-(1-delta))/alpha)^(1/(alpha-1));
+c = k^alpha-delta*k;
+z = 0;
+U = ln(c);
+W = U/(1-beta);
+end;
+
+shocks;
+var e;
+stderr 1;
+end;
+
+steady;
+resid;
+stoch_simul(order=2, irf=0);
diff --git a/tests/optimal_policy/neo_growth_ramsey.mod b/tests/optimal_policy/neo_growth_ramsey.mod
new file mode 100644
index 0000000000..aa2e637d6c
--- /dev/null
+++ b/tests/optimal_policy/neo_growth_ramsey.mod
@@ -0,0 +1,85 @@
+/*
+ * 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 <http://www.gnu.org/licenses/>.
+ */
+/*
+ * This file computes a second-order approximation of the neo-classical growth model.
+ * It assesses the conditional and unconditional welfares computed by the evaluate_planner_objective function
+ * and compares them to a by-hand assessment stemming from the results the model neo_growth.mod incur.
+ */
+
+var k z c;
+
+varexo e;
+
+parameters beta gamma alpha delta rho s;
+
+beta = 0.987;
+gamma = 1;
+delta = 0.012;
+alpha = 0.4;
+rho = 0.95;
+s = 0.007;
+
+model;
+k=exp(z)*k(-1)^(alpha)-c+(1-delta)*k(-1);
+z=rho*z(-1)+s*e;
+end;
+
+steady_state_model;
+z = 0;
+end;
+
+shocks;
+var e;
+stderr 1;
+end;
+
+planner_objective ln(c);
+ramsey_model(instruments=(k,c), planner_discount=beta);
+
+initval;
+   k = ((1/beta-(1-delta))/alpha)^(1/(alpha-1));
+   c = k^alpha-delta*k;
+end;
+
+steady;
+resid;
+stoch_simul(order=2, irf=0);
+
+planner_objective_value = evaluate_planner_objective(M_, options_, oo_);
+
+if ~exist('neo_growth_results.mat','file');
+   error('neo_growth must be run first');
+end;
+
+oo1 = load('neo_growth_results','oo_');
+M1 = load('neo_growth_results','M_');
+options1 = load('neo_growth_results','options_');
+unc_W_hand = oo1.oo_.mean(strmatch('W',M1.M_.endo_names,'exact'));
+
+initial_condition_states = repmat(oo1.oo_.dr.ys,1,M1.M_.maximum_lag);
+shock_matrix = zeros(1,M1.M_.exo_nbr);
+y_sim = simult_(M1.M_,options1.options_,initial_condition_states,oo1.oo_.dr,shock_matrix,options1.options_.order);
+cond_W_hand=y_sim(strmatch('W',M1.M_.endo_names,'exact'),2);
+
+if abs((unc_W_hand - planner_objective_value(1))/unc_W_hand) > 1e-6;
+   error('Inaccurate unconditional welfare assessment');
+end;
+if abs((cond_W_hand - planner_objective_value(2))/cond_W_hand) > 1e-6;
+   error('Inaccurate conditional welfare assessment');
+end;
-- 
GitLab