From 6f1c02769e743d9f3c81b824dce7c521861d7db9 Mon Sep 17 00:00:00 2001 From: Michel Juillard Date: Tue, 24 May 2011 15:33:09 +0200 Subject: [PATCH] fixing evaluation of objective function under optimal policy --- matlab/evaluate_planner_objective.m | 126 +++++++++++++++------------- matlab/ramsey_policy.m | 2 +- 2 files changed, 70 insertions(+), 58 deletions(-) diff --git a/matlab/evaluate_planner_objective.m b/matlab/evaluate_planner_objective.m index cda2169a5..a393ce1ec 100644 --- a/matlab/evaluate_planner_objective.m +++ b/matlab/evaluate_planner_objective.m @@ -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 . -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 \ No newline at end of file diff --git a/matlab/ramsey_policy.m b/matlab/ramsey_policy.m index 259d6787d..ea778d2d6 100644 --- a/matlab/ramsey_policy.m +++ b/matlab/ramsey_policy.m @@ -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; \ No newline at end of file