Merge branch 'perfect_foresight_planner' into 'master'

Ref. #1680: perfect foresight

See merge request Dynare/dynare!1833

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;
-        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;
+    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
-
-        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

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 \

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);

tests/optimal_policy/neo_growth_common.inc

+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;

tests/optimal_policy/neo_growth_foresight.mod

+/*
+ * 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 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

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_);

tests/optimal_policy/neo_growth_ramsey_common.inc

+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;

tests/optimal_policy/neo_growth_ramsey_foresight.mod

+/*
+ * 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 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;