fixing evaluation of objective function under optimal policy
Michel Juillard
parent
ae56c42b96
commit
6f1c02769e
|
|
@ -1,4 +1,4 @@
|
|||
function oo1 = evaluate_planner_objective(dr,M,oo,options)
|
||||
function planner_objective_value = evaluate_planner_objective(M,oo,options)
|
||||
|
||||
%function oo1 = evaluate_planner_objective(dr,M,oo,options)
|
||||
% computes value of planner objective function
|
||||
|
|
@ -32,87 +32,99 @@ function oo1 = evaluate_planner_objective(dr,M,oo,options)
|
|||
% You should have received a copy of the GNU General Public License
|
||||
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
oo1 = oo;
|
||||
dr = oo.dr;
|
||||
endo_nbr = M.endo_nbr;
|
||||
exo_nbr = M.exo_nbr;
|
||||
nstatic = dr.nstatic;
|
||||
npred = dr.npred;
|
||||
lead_lag_incidence = M.lead_lag_incidence;
|
||||
if (size(lead_lag_incidence,1) == 5) ...
|
||||
&& (nnz(lead_lag_incidence(1,:)) > 0 ...
|
||||
|| nnz(lead_lag_incidence(5,:)) > 0)
|
||||
error(['procedure not yet implemented for leads and lags on more ' ...
|
||||
'than one period'])
|
||||
end
|
||||
beta = get_optimal_policy_discount_factor(M.params,M.param_names);
|
||||
if options.ramsey_policy
|
||||
i_org = (1:M.orig_model.endo_nbr)';
|
||||
i_org = (1:M.orig_endo_nbr)';
|
||||
else
|
||||
i_org = (1:M.endo_nbr)';
|
||||
end
|
||||
ipred = find(lead_lag_incidence(M.maximum_lag,:))';
|
||||
orig_ifrwrd = find(M.orig_model.lead_lag_incidence(3,:));
|
||||
orig_npred = nnz(M.orig_model.lead_lag_incidence(1,:));
|
||||
order_var = dr.order_var;
|
||||
ys = oo.dr.ys;
|
||||
[U,Uy,Uyy] = feval([M.fname '_objective_static'],ys,zeros(1,exo_nbr), ...
|
||||
M.params);
|
||||
z = repmat(ys(1:M.orig_model.endo_nbr),1,3);
|
||||
z = z(find(M.orig_model.lead_lag_incidence'));
|
||||
[F,Fy,Fyy] = feval([M.fname '_dynamic'],z,zeros(3,exo_nbr), ...
|
||||
M.params,2);
|
||||
mu_ss = oo.dr.ys(M.orig_model.endo_nbr+exo_nbr+(1:size(F,1)));
|
||||
%K = reshape(1:endo_nbr^2,endo_nbr,endo_nbr);
|
||||
%K = K(order_var,order_var);
|
||||
%Uyy = Uyy(:,K(:));
|
||||
beta = options.planner_discount;
|
||||
LQ = true;
|
||||
|
||||
Gy = dr.ghx(nstatic+(1:npred),:);
|
||||
Gu = dr.ghu(nstatic+(1:npred),:);
|
||||
gy(dr.order_var,:) = dr.ghx;
|
||||
gu(dr.order_var,:) = dr.ghu;
|
||||
gy = gy(i_org,:);
|
||||
gu = gu(i_org,:);
|
||||
muFy = mu_ss'*Fy;
|
||||
muFyy = mu_ss'*Fyy;
|
||||
Zy = [eye(orig_npred,npred); gy; gy(orig_ifrwrd,:)*Gy; zeros(exo_nbr, ...
|
||||
npred)];
|
||||
Zu = [zeros(orig_npred,exo_nbr); gu; gy(orig_ifrwrd,:)*Gu; eye(exo_nbr, ...
|
||||
exo_nbr)];
|
||||
Zyy = kron(Zy,Zy);
|
||||
Zuu = kron(Zu,Zu);
|
||||
Zyu = kron(Zy,Zu);
|
||||
|
||||
Wbar = U/(1-beta);
|
||||
Wy = (Uy*gy+muFy*Zy)/(eye(npred)-beta*Gy);
|
||||
Wu = Uy*gu + muFy*Zu + beta*Wy*Gu;
|
||||
% N1 = Uyy*kron(gy,gy)
|
||||
% N2 = muFyy*Zyy
|
||||
format long
|
||||
Wyy = (Uyy*kron(gy,gy)+muFyy*Zyy)/(eye(npred^2)-beta*kron(Gy,Gy))
|
||||
Wyu = Uyy*kron(gy,gu) + muFyy*Zyu + beta*Wyy*kron(Gy,Gu);
|
||||
Wuu = Uyy*kron(gu,gu) + muFyy*Zuu + beta*Wyy*kron(Gu,Gu);
|
||||
Wss = beta*Wuu*vec(M.Sigma_e)/(1-beta);
|
||||
if options.ramsey_policy && options.order == 1 && ~options.linear
|
||||
options.order = 2;
|
||||
options.qz_criterium = 1+1e-6;
|
||||
[dr,info] = dr1(oo.dr,0,M,options,oo);
|
||||
Gyy = dr.ghxx(nstatic+(1:npred),:);
|
||||
Guu = dr.ghuu(nstatic+(1:npred),:);
|
||||
Gyu = dr.ghxu(nstatic+(1:npred),:);
|
||||
Gss = dr.ghs2(nstatic+(1:npred),:);
|
||||
gyy(dr.order_var,:) = dr.ghxx;
|
||||
guu(dr.order_var,:) = dr.ghuu;
|
||||
gyu(dr.order_var,:) = dr.ghxu;
|
||||
gss(dr.order_var,:) = dr.ghs2;
|
||||
LQ = false;
|
||||
end
|
||||
|
||||
ys = oo.dr.ys;
|
||||
|
||||
u = oo.exo_simul(1,:)';
|
||||
|
||||
[U,Uy,Uyy] = feval([M.fname '_objective_static'],ys,zeros(1,exo_nbr), ...
|
||||
M.params);
|
||||
|
||||
Uyy = full(Uyy);
|
||||
|
||||
[err,Uyygygy] = A_times_B_kronecker_C(Uyy,gy,gy,options.threads.kronecker.A_times_B_kronecker_C);
|
||||
mexErrCheck('A_times_B_kronecker_C', err);
|
||||
[err,Uyygugu] = A_times_B_kronecker_C(Uyy,gu,gu,options.threads.kronecker.A_times_B_kronecker_C);
|
||||
mexErrCheck('A_times_B_kronecker_C', err);
|
||||
[err,Uyygygu] = A_times_B_kronecker_C(Uyy,gy,gu,options.threads.kronecker.A_times_B_kronecker_C);
|
||||
mexErrCheck('A_times_B_kronecker_C', err);
|
||||
|
||||
Wbar =U/(1-beta);
|
||||
Wy = Uy*gy/(eye(npred)-beta*Gy);
|
||||
Wu = Uy*gu+beta*Wy*Gu;
|
||||
if LQ
|
||||
Wyy = Uyygygy/(eye(npred*npred)-beta*kron(Gy,Gy));
|
||||
else
|
||||
Wyy = (Uy*gyy+Uyygygy+beta*Wy*Gyy)/(eye(npred*npred)-beta*kron(Gy,Gy));
|
||||
end
|
||||
[err,Wyygugu] = A_times_B_kronecker_C(Wyy,Gu,Gu,options.threads.kronecker.A_times_B_kronecker_C);
|
||||
mexErrCheck('A_times_B_kronecker_C', err);
|
||||
[err,Wyygygu] = A_times_B_kronecker_C(Wyy,Gy,Gu,options.threads.kronecker.A_times_B_kronecker_C);
|
||||
mexErrCheck('A_times_B_kronecker_C', err);
|
||||
if LQ
|
||||
Wuu = Uyygugu+beta*Wyygugu;
|
||||
Wyu = Uyygygu+beta*Wyygygu;
|
||||
Wss = beta*Wuu*M.Sigma_e(:)/(1-beta);
|
||||
else
|
||||
Wuu = Uy*guu+Uyygugu+beta*(Wy*Guu+Wyygugu);
|
||||
Wyu = Uy*gyu+Uyygygu+beta*(Wy*Gyu+Wyygygu);
|
||||
Wss = (Uy*gss+beta*(Wuu*M.Sigma_e(:)+Wy*Gss))/(1-beta);
|
||||
end
|
||||
if options.ramsey_policy
|
||||
yhat = [oo.endo_simul; ...
|
||||
zeros(M.endo_nbr-size(oo.endo_simul,1),1)];
|
||||
else
|
||||
yhat = oo.endo_simul;
|
||||
end
|
||||
yhat = yhat(dr.order_var(nstatic+(1:npred)),1)-dr.ys(dr.order_var(nstatic+(1:npred)))
|
||||
yhat = yhat(dr.order_var(nstatic+(1:npred)),1)-dr.ys(dr.order_var(nstatic+(1:npred)));
|
||||
u = oo.exo_simul(1,:)';
|
||||
|
||||
planner_objective_value = Wbar+Wy*yhat+Wu*u+Wyu*kron(yhat,u) ...
|
||||
+ 0.5*(Wyy*kron(yhat,yhat) + Wuu*kron(u,u)+Wss);
|
||||
disp(Wyy*kron(yhat,yhat))
|
||||
oo1.planner_objective_value = planner_objective_value;
|
||||
[err,Wyyyhatyhat] = A_times_B_kronecker_C(Wyy,yhat,yhat,options.threads.kronecker.A_times_B_kronecker_C);
|
||||
mexErrCheck('A_times_B_kronecker_C', err);
|
||||
[err,Wuuuu] = A_times_B_kronecker_C(Wuu,u,u,options.threads.kronecker.A_times_B_kronecker_C);
|
||||
mexErrCheck('A_times_B_kronecker_C', err);
|
||||
[err,Wyuyhatu] = A_times_B_kronecker_C(Wyu,yhat,u,options.threads.kronecker.A_times_B_kronecker_C);
|
||||
mexErrCheck('A_times_B_kronecker_C', err);
|
||||
planner_objective_value = Wbar+Wy*yhat+Wu*u+Wyuyhatu ...
|
||||
+ 0.5*(Wyyyhatyhat + Wuuuu+Wss);
|
||||
if ~options.noprint
|
||||
disp(' ')
|
||||
if options.ramsey_policy
|
||||
disp(sprintf('Value of planner objective function under Ramsey policy: %16.8f', ...
|
||||
planner_objective_value))
|
||||
else
|
||||
disp(['Value of planner objective function: ' ...
|
||||
num2str(planner_objective_value)])
|
||||
disp(' ')
|
||||
end
|
||||
disp(['Approximated value of planner objective function: ' ...
|
||||
num2str(planner_objective_value)])
|
||||
disp(' ')
|
||||
end
|
||||
|
|
@ -35,6 +35,6 @@ if options_.noprint == 0
|
|||
end
|
||||
|
||||
|
||||
%oo_ = evaluate_planner_objective(oo_.dr,M_,oo_,options_);
|
||||
oo_.planner_objective_value = evaluate_planner_objective(M_,oo_,options_);
|
||||
|
||||
options_ = oldoptions;
|
||||
Loading…
Reference in New Issue