evaluate_planner_objective

function oo1 = evaluate_planner_objective(dr,M,oo,options)
  computes value of planner objective function     
 
 INPUTS
   dr:       (structure) decision rule
   M:        (structure) model description
   oo:       (structure) output results
   options:  (structure) options
    
 OUTPUTS
   oo1:      (structure) updated output results

 SPECIAL REQUIREMENTS
   none

0001 function planner_objective_value = evaluate_planner_objective(M,oo,options)
0002 
0003 %function oo1 = evaluate_planner_objective(dr,M,oo,options)
0004 %  computes value of planner objective function
0005 %
0006 % INPUTS
0007 %   dr:       (structure) decision rule
0008 %   M:        (structure) model description
0009 %   oo:       (structure) output results
0010 %   options:  (structure) options
0011 %
0012 % OUTPUTS
0013 %   oo1:      (structure) updated output results
0014 %
0015 % SPECIAL REQUIREMENTS
0016 %   none
0017 
0018 % Copyright (C) 2007-2011 Dynare Team
0019 %
0020 % This file is part of Dynare.
0021 %
0022 % Dynare is free software: you can redistribute it and/or modify
0023 % it under the terms of the GNU General Public License as published by
0024 % the Free Software Foundation, either version 3 of the License, or
0025 % (at your option) any later version.
0026 %
0027 % Dynare is distributed in the hope that it will be useful,
0028 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0029 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0030 % GNU General Public License for more details.
0031 %
0032 % You should have received a copy of the GNU General Public License
0033 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0034 
0035 dr = oo.dr;
0036 endo_nbr = M.endo_nbr;
0037 exo_nbr = M.exo_nbr;
0038 nstatic = dr.nstatic;
0039 npred = dr.npred;
0040 lead_lag_incidence = M.lead_lag_incidence;
0041 beta = get_optimal_policy_discount_factor(M.params,M.param_names);
0042 if options.ramsey_policy
0043     i_org = (1:M.orig_endo_nbr)';
0044 else
0045     i_org = (1:M.endo_nbr)';
0046 end
0047 ipred = find(lead_lag_incidence(M.maximum_lag,:))';
0048 order_var = dr.order_var;
0049 LQ = true;
0050 
0051 Gy = dr.ghx(nstatic+(1:npred),:);
0052 Gu = dr.ghu(nstatic+(1:npred),:);
0053 gy(dr.order_var,:) = dr.ghx;
0054 gu(dr.order_var,:) = dr.ghu;
0055 
0056 if options.ramsey_policy && options.order == 1 && ~options.linear
0057     options.order = 2;
0058     options.qz_criterium = 1+1e-6;
0059     [dr,info] = stochastic_solvers(oo.dr,0,M,options,oo);
0060     Gyy = dr.ghxx(nstatic+(1:npred),:);
0061     Guu = dr.ghuu(nstatic+(1:npred),:);
0062     Gyu = dr.ghxu(nstatic+(1:npred),:);
0063     Gss = dr.ghs2(nstatic+(1:npred),:);
0064     gyy(dr.order_var,:) = dr.ghxx;
0065     guu(dr.order_var,:) = dr.ghuu;
0066     gyu(dr.order_var,:) = dr.ghxu;
0067     gss(dr.order_var,:) = dr.ghs2;
0068     LQ = false;
0069 end
0070 
0071 ys = oo.dr.ys;
0072 
0073 u = oo.exo_simul(1,:)';
0074 
0075 [U,Uy,Uyy] = feval([M.fname '_objective_static'],ys,zeros(1,exo_nbr), ...
0076                    M.params);
0077 
0078 Uyy = full(Uyy);
0079 
0080 [Uyygygy, err] = A_times_B_kronecker_C(Uyy,gy,gy,options.threads.kronecker.A_times_B_kronecker_C);
0081 mexErrCheck('A_times_B_kronecker_C', err);
0082 [Uyygugu, err] = A_times_B_kronecker_C(Uyy,gu,gu,options.threads.kronecker.A_times_B_kronecker_C);
0083 mexErrCheck('A_times_B_kronecker_C', err);
0084 [Uyygygu, err] = A_times_B_kronecker_C(Uyy,gy,gu,options.threads.kronecker.A_times_B_kronecker_C);
0085 mexErrCheck('A_times_B_kronecker_C', err);
0086 
0087 Wbar =U/(1-beta);
0088 Wy = Uy*gy/(eye(npred)-beta*Gy);
0089 Wu = Uy*gu+beta*Wy*Gu;
0090 if LQ
0091     Wyy = Uyygygy/(eye(npred*npred)-beta*kron(Gy,Gy));
0092 else
0093     Wyy = (Uy*gyy+Uyygygy+beta*Wy*Gyy)/(eye(npred*npred)-beta*kron(Gy,Gy));
0094 end
0095 [Wyygugu, err] = A_times_B_kronecker_C(Wyy,Gu,Gu,options.threads.kronecker.A_times_B_kronecker_C);
0096 mexErrCheck('A_times_B_kronecker_C', err);
0097 [Wyygygu,err] = A_times_B_kronecker_C(Wyy,Gy,Gu,options.threads.kronecker.A_times_B_kronecker_C);
0098 mexErrCheck('A_times_B_kronecker_C', err);
0099 if LQ
0100     Wuu = Uyygugu+beta*Wyygugu;
0101     Wyu = Uyygygu+beta*Wyygygu;
0102     Wss = beta*Wuu*M.Sigma_e(:)/(1-beta);
0103 else
0104     Wuu = Uy*guu+Uyygugu+beta*(Wy*Guu+Wyygugu);
0105     Wyu = Uy*gyu+Uyygygu+beta*(Wy*Gyu+Wyygygu);
0106     Wss = (Uy*gss+beta*(Wuu*M.Sigma_e(:)+Wy*Gss))/(1-beta);
0107 end
0108 if options.ramsey_policy
0109     yhat = zeros(M.endo_nbr,1);
0110     yhat(1:M.orig_endo_nbr) = oo.steady_state(1:M.orig_endo_nbr);
0111 else
0112     yhat = oo.endo_simul;
0113 end
0114 yhat = yhat(dr.order_var(nstatic+(1:npred)),1)-dr.ys(dr.order_var(nstatic+(1:npred)));
0115 u = oo.exo_simul(1,:)';
0116 
0117 [Wyyyhatyhat, err] = A_times_B_kronecker_C(Wyy,yhat,yhat,options.threads.kronecker.A_times_B_kronecker_C);
0118 mexErrCheck('A_times_B_kronecker_C', err);
0119 [Wuuuu, err] = A_times_B_kronecker_C(Wuu,u,u,options.threads.kronecker.A_times_B_kronecker_C);
0120 mexErrCheck('A_times_B_kronecker_C', err);
0121 [Wyuyhatu, err] = A_times_B_kronecker_C(Wyu,yhat,u,options.threads.kronecker.A_times_B_kronecker_C);
0122 mexErrCheck('A_times_B_kronecker_C', err);
0123 planner_objective_value(1) = Wbar+Wy*yhat+Wu*u+Wyuyhatu ...
0124     + 0.5*(Wyyyhatyhat + Wuuuu+Wss);
0125 planner_objective_value(2) = Wbar + 0.5*Wss;
0126 if ~options.noprint
0127     disp(' ')
0128     disp('Approximated value of planner objective function')
0129     disp(['    - with initial Lagrange multipliers set to 0: ' ...
0130           num2str(planner_objective_value(1)) ])
0131     disp(['    - with initial Lagrange multipliers set to steady state: ' ...
0132           num2str(planner_objective_value(2)) ])
0133     disp(' ')
0134 end