diff --git a/matlab/evaluate_planner_objective.m b/matlab/evaluate_planner_objective.m
index 878b06ec3..3f09ca3fa 100644
--- a/matlab/evaluate_planner_objective.m
+++ b/matlab/evaluate_planner_objective.m
@@ -49,6 +49,8 @@ function planner_objective_value = evaluate_planner_objective(M_,options_,oo_)
% E(W) = (1 - beta)^{-1} ( Ubar + 0.5 ( U_xx h_y^2 E(yhat^2) + U_xx h_u^2 E(u^2) )
% As for the conditional welfare, the second-order formula above is still valid, but the derivatives of W no longer contain any second-order derivatives of the policy and transition functions h and g.
+% In the deterministic case, resorting to approximations for welfare is no longer required as it is possible to simulate the model given initial conditions for pre-determined variables and terminal conditions for forward-looking variables, whether these initial and terminal conditions are explicitly or implicitly specified. Assuming that the number of simulated periods is high enough for the new steady-state to be reached, the new unconditional welfare is thus the last period's welfare. As for the conditional welfare, it can be derived using backward recursions on the equation W = U + beta*W(+1) starting from the final unconditional steady-state welfare.
+
% INPUTS
% M_: (structure) model description
% options_: (structure) options
@@ -81,87 +83,101 @@ 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_.ramsey_policy
- 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;
+ if oo_.gui.ran_perfect_foresight
+ T = size(oo_.endo_simul,2);
+ [U_term] = feval([M_.fname '.objective.static'],oo_.endo_simul(:,T),oo_.exo_simul(T,:), M_.params);
+ EW = U_term/(1-beta);
+ W = EW;
+ for t=T:-1:2
+ [U] = feval([M_.fname '.objective.static'],oo_.endo_simul(:,t),oo_.exo_simul(t,:), M_.params);
+ W = U + beta*W;
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)));
-
- 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
+ ys = oo_.dr.ys;
+ 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
+ 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)));
+
+ 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
end
elseif options_.discretionary_policy
-
+ ys = oo_.dr.ys;
+
[U,Uy,Uyy] = feval([M_.fname '.objective.static'],ys,zeros(1,exo_nbr), M_.params);
Gy = dr.ghx(nstatic+(1:nspred),:);
@@ -221,11 +237,15 @@ elseif options_.discretionary_policy
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))
+ if oo_.gui.ran_perfect_foresight
+ fprintf('\nSimulated value of unconditional welfare: %10.8f\n', planner_objective_value(1))
+ fprintf('\nSimulated value of conditional welfare: %10.8f\n', planner_objective_value(2))
+ else
+ 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))
+ end
elseif options_.discretionary_policy
fprintf('\nApproximated value of unconditional welfare with discretionary policy: %10.8f\n', planner_objective_value(1))
fprintf('\nApproximated value of conditional welfare with discretionary policy: %10.8f\n', planner_objective_value(2))
-
end
end
diff --git a/tests/Makefile.am b/tests/Makefile.am
index a0a36a1d4..14e6c86f2 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -580,6 +580,9 @@ MFILES = histval_initval_file/ramst_initval_file_data.m
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
+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
+
example1_use_dll.m.trs: example1.m.trs
example1_use_dll.o.trs: example1.o.trs
@@ -1185,6 +1188,8 @@ EXTRA_DIST = \
observation_trends_and_prefiltering/Trend_model_calib_no_prefilter_common.inc \
observation_trends_and_prefiltering/Trend_load_data_common.inc \
observation_trends_and_prefiltering/Trend_no_prefilter_conditional_forecast.inc \
+ optimal_policy/neo_growth_common.inc \
+ optimal_policy/neo_growth_ramsey_common.inc \
optimal_policy/Ramsey/oo_ramsey_policy_initval.mat \
optimizers/optimizer_function_wrapper.m \
optimizers/fs2000.common.inc \
diff --git a/tests/optimal_policy/neo_growth.mod b/tests/optimal_policy/neo_growth.mod
index 2c5165c8c..2b6afa29e 100644
--- a/tests/optimal_policy/neo_growth.mod
+++ b/tests/optimal_policy/neo_growth.mod
@@ -21,34 +21,7 @@
* 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;
+@#include "neo_growth_common.inc"
shocks;
var e;
@@ -57,4 +30,5 @@ end;
steady;
resid;
+
stoch_simul(order=2, irf=0);
diff --git a/tests/optimal_policy/neo_growth_common.inc b/tests/optimal_policy/neo_growth_common.inc
new file mode 100644
index 000000000..0d5d1d6b3
--- /dev/null
+++ b/tests/optimal_policy/neo_growth_common.inc
@@ -0,0 +1,28 @@
+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;
diff --git a/tests/optimal_policy/neo_growth_foresight.mod b/tests/optimal_policy/neo_growth_foresight.mod
new file mode 100644
index 000000000..c1b43937f
--- /dev/null
+++ b/tests/optimal_policy/neo_growth_foresight.mod
@@ -0,0 +1,41 @@
+/*
+ * 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 .
+ */
+ /*
+ * This file simulates the neo-classical growth model in a perfect foresight framework.
+ * It is called by neo_growth_ramsey_foresight.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"
+
+initval;
+z = 0;
+end;
+
+steady;
+
+shocks;
+var e;
+periods 1;
+values 1;
+end;
+
+resid;
+
+perfect_foresight_setup(periods=200);
+perfect_foresight_solver;
\ No newline at end of file
diff --git a/tests/optimal_policy/neo_growth_ramsey.mod b/tests/optimal_policy/neo_growth_ramsey.mod
index aa2e637d6..2258e7012 100644
--- a/tests/optimal_policy/neo_growth_ramsey.mod
+++ b/tests/optimal_policy/neo_growth_ramsey.mod
@@ -22,43 +22,13 @@
* 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;
+@#include "neo_growth_ramsey_common.inc"
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_);
diff --git a/tests/optimal_policy/neo_growth_ramsey_common.inc b/tests/optimal_policy/neo_growth_ramsey_common.inc
new file mode 100644
index 000000000..01e4a6d39
--- /dev/null
+++ b/tests/optimal_policy/neo_growth_ramsey_common.inc
@@ -0,0 +1,32 @@
+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;
+
+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;
diff --git a/tests/optimal_policy/neo_growth_ramsey_foresight.mod b/tests/optimal_policy/neo_growth_ramsey_foresight.mod
new file mode 100644
index 000000000..df9c5ef2a
--- /dev/null
+++ b/tests/optimal_policy/neo_growth_ramsey_foresight.mod
@@ -0,0 +1,53 @@
+/*
+ * 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 .
+ */
+ /*
+ * This file simulates a perfect-foresight version 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 of the model neo_growth_foresight.mod
+ */
+
+@#include "neo_growth_ramsey_common.inc"
+
+shocks;
+var e;
+periods 1;
+values 1;
+end;
+
+perfect_foresight_setup(periods=200);
+perfect_foresight_solver;
+
+planner_objective_value = evaluate_planner_objective(M_, options_, oo_);
+
+if ~exist('neo_growth_foresight_results.mat','file');
+ error('neo_growth_foresight must be run first');
+end;
+
+oo1 = load('neo_growth_foresight_results','oo_');
+M1 = load('neo_growth_foresight_results','M_');
+options1 = load('neo_growth_foresight_results','options_');
+cond_W_hand = oo1.oo_.endo_simul(strmatch('W',M1.M_.endo_names,'exact'),2);
+unc_W_hand = oo1.oo_.endo_simul(strmatch('W',M1.M_.endo_names,'exact'),end);
+
+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;