Skip to content
GitLab
Commit e880d1bc authored by Normann Rion's avatar Normann Rion Committed by Normann Rion
Browse files

Ref. #1680: 2nd-order welfare

parent 5485597a
Hide whitespace changes
Inline Side-by-side
matlab/evaluate_planner_objective.m
View file @ e880d1bc
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
matlab/ramsey_policy.m
View file @ e880d1bc
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);
......
tests/Makefile.am
View file @ e880d1bc
......@@ -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
......
tests/optimal_policy/Ramsey/ramsey_ex.mod
View file @ e880d1bc
......@@ -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
tests/optimal_policy/Ramsey/ramsey_ex_initval.mod
View file @ e880d1bc
......@@ -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 ...
......
tests/optimal_policy/neo_growth.mod 0 → 100644
View file @ e880d1bc
/*
* 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);
tests/optimal_policy/neo_growth_ramsey.mod 0 → 100644
View file @ e880d1bc
/*
* 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;
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment