diff --git a/matlab/evaluate_planner_objective.m b/matlab/evaluate_planner_objective.m
index a6df353260b7916ce77032f1a111f906a366be39..26e23b4a15e1533faa7ae35e8a7aee99f0e9517b 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 720e2d83374c929e347f3a55d7c37f57b06430c3..69b2a03b608764ccc773b803d2bcd7e600017fb8 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 46f29f536763a3a01d24703dc1f7e063a6870a13..07dcebb9b63dcef2b1a9ff335ea2562cf3f74459 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
Binary files a/tests/optimal_policy/Ramsey/oo_ramsey_policy_initval.mat and b/tests/optimal_policy/Ramsey/oo_ramsey_policy_initval.mat differ
diff --git a/tests/optimal_policy/Ramsey/ramsey_ex.mod b/tests/optimal_policy/Ramsey/ramsey_ex.mod
index d28882f970334f18e752b5c89b1110022200438d..04e124462638f8527c70184ad60e40ee6ba1fea6 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 170b956f99aca58e2d75bda91de8de2e973598ff..9f0d8e1326b6a7f28ad949345b9b44d3e510ae27 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 0000000000000000000000000000000000000000..2c5165c8c7706a20dc504d02c626c195d2b31bd7
--- /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 0000000000000000000000000000000000000000..aa2e637d6cc651eeee76f7a93ea7edc9a50cebd4
--- /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;