From cc3aeafd00f5135425dcc869e54a9c99c1e20790 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Wed, 26 Aug 2015 10:00:35 +0200 Subject: [PATCH 1/7] Save planner objective function with discretionary_policy --- matlab/discretionary_policy.m | 2 +- matlab/evaluate_planner_objective.m | 12 ++++++++---- 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/matlab/discretionary_policy.m b/matlab/discretionary_policy.m index eda6aa11d..d12af30bd 100644 --- a/matlab/discretionary_policy.m +++ b/matlab/discretionary_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_,options_,oo_); options_ = oldoptions; \ No newline at end of file diff --git a/matlab/evaluate_planner_objective.m b/matlab/evaluate_planner_objective.m index 32db62625..e6f80f9bd 100644 --- a/matlab/evaluate_planner_objective.m +++ b/matlab/evaluate_planner_objective.m @@ -92,7 +92,6 @@ if ~isempty(M.endo_histval) end end yhat1 = yhat1(dr.order_var(nstatic+(1:nspred)),1)-dr.ys(dr.order_var(nstatic+(1:nspred))); -yhat2 = yhat2(dr.order_var(nstatic+(1:nspred)),1)-dr.ys(dr.order_var(nstatic+(1:nspred))); u = oo.exo_simul(1,:)'; [Wyyyhatyhat1, err] = A_times_B_kronecker_C(Wyy,yhat1,yhat1,options.threads.kronecker.A_times_B_kronecker_C); @@ -104,6 +103,7 @@ mexErrCheck('A_times_B_kronecker_C', err); 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, err] = A_times_B_kronecker_C(Wyy,yhat2,yhat2,options.threads.kronecker.A_times_B_kronecker_C); mexErrCheck('A_times_B_kronecker_C', err); [Wyuyhatu2, err] = A_times_B_kronecker_C(Wyu,yhat2,u,options.threads.kronecker.A_times_B_kronecker_C); @@ -115,9 +115,13 @@ end if ~options.noprint skipline() disp('Approximated value of planner objective function') - disp([' - with initial Lagrange multipliers set to 0: ' ... + if options.ramsey_policy + disp([' - with initial Lagrange multipliers set to 0: ' ... num2str(planner_objective_value(2)) ]) - disp([' - with initial Lagrange multipliers set to steady state: ' ... - num2str(planner_objective_value(1)) ]) + disp([' - with initial Lagrange multipliers set to steady state: ' ... + num2str(planner_objective_value(1)) ]) + elseif options.discretionary_policy + fprintf('with discretionary policy: %10.8f',planner_objective_value(1)) + end skipline() end From c66ee8941ae8b81f3e0b177fbf7ba9fe1d121ed6 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Sun, 30 Aug 2015 11:14:35 +0200 Subject: [PATCH 2/7] Document discretionary_policy_engine.m --- matlab/discretionary_policy_engine.m | 72 ++++++++++++++++++++++------ 1 file changed, 57 insertions(+), 15 deletions(-) diff --git a/matlab/discretionary_policy_engine.m b/matlab/discretionary_policy_engine.m index 9b94d7301..154b66426 100644 --- a/matlab/discretionary_policy_engine.m +++ b/matlab/discretionary_policy_engine.m @@ -1,16 +1,54 @@ function [H,G,retcode]=discretionary_policy_engine(AAlag,AA0,AAlead,BB,bigw,instr_id,beta,solve_maxit,discretion_tol,qz_criterium,H00,verbose) % Solves the discretionary problem for a model of the form: -% AAlag*yy_{t-1}+AA0*yy_t+AAlead*yy_{t+1}+BB*e=0, with W the weight on the -% variables in vector y_t and instr_id is the location of the instruments -% in the yy_t vector. +% +% Loss=E_0 sum_{t=0}^{\infty} beta^t [y_t'*W*y+x_t'*Q*x_t] +% subject to +% AAlag*yy_{t-1}+AA0*yy_t+AAlead*yy_{t+1}+BB*e=0 +% +% with W the weight on the variables in vector y_t. +% +% The solution takes the form +% y_t=H*y_{t-1}+G*e_t +% where H=[H1;F1] and G=[H2;F2]. +% % We use the Dennis (2007, Macroeconomic Dynamics) algorithm and so we need % to re-write the model in the form % A0*y_t=A1*y_{t-1}+A2*y_{t+1}+A3*x_t+A4*x_{t+1}+A5*e_t, with W the % weight on the y_t vector and Q the weight on the x_t vector of % instruments. +% +% Inputs: +% AAlag [double] matrix of coefficients on lagged +% variables +% AA0 [double] matrix of coefficients on +% contemporaneous variables +% AAlead [double] matrix of coefficients on +% leaded variables +% BB [double] matrix of coefficients on +% shocks +% bigw [double] matrix of coefficients on variables in +% loss/objective function; stacks [W and Q] +% instr_id [double] location vector of the instruments in the yy_t vector. +% beta [scalar] planner discount factor +% solve_maxit [scalar] maximum number of iterations +% discretion_tol [scalar] convergence criterion for solution +% qz_criterium [scalar] tolerance for QZ decomposition +% H00 +% verbose [scalar] dummy to control verbosity +% +% Outputs: +% H [double] (endo_nbr*endo_nbr) solution matrix for endogenous +% variables, stacks [H1 and H1] +% G [double] (endo_nbr*exo_nbr) solution matrix for shocks, stacks [H2 and F2] +% +% retcode [scalar] return code +% +% Algorithm: +% Dennis, Richard (2007): Optimal policy in rational expectations models: new solution algorithms, +% Macroeconomic Dynamics, 11, 31–55. -% Copyright (C) 2007-2012 Dynare Team +% Copyright (C) 2007-2015 Dynare Team % % This file is part of Dynare. % @@ -53,14 +91,15 @@ end [A0,A1,A2,A3,A4,A5,W,Q,endo_nbr,exo_nbr,aux,endo_augm_id]=GetDennisMatrices(AAlag,AA0,AAlead,BB,bigw,instr_id); % aux is a logical index of the instruments which appear with lags in the -% model. Their location in the state vector is instr_id(aux) +% model. Their location in the state vector is instr_id(aux); % endo_augm_id is index (not logical) of locations of the augmented vector % of non-instrumental variables AuxiliaryVariables_nbr=sum(aux); H0=zeros(endo_nbr+AuxiliaryVariables_nbr); if ~isempty(H00) - H0(1:endo_nbr,1:endo_nbr)=H00;clear H00 + H0(1:endo_nbr,1:endo_nbr)=H00; + clear H00 end H10=H0(endo_augm_id,endo_augm_id); @@ -69,6 +108,7 @@ F10=H0(instr_id,endo_augm_id); iter=0; H1=H10; F1=F10; +%solve equations (20) and (22) via fixed point iteration while 1 iter=iter+1; P=SylvesterDoubling(W+beta*F1'*Q*F1,beta*H1',H1,discretion_tol,solve_maxit); @@ -79,11 +119,11 @@ while 1 return end end - D=A0-A2*H1-A4*F1; + D=A0-A2*H1-A4*F1; %equation (20) Dinv=inv(D); - A3DPD=A3'*Dinv'*P*Dinv; - F1=-(Q+A3DPD*A3)\(A3DPD*A1); - H1=Dinv*(A1+A3*F1); + A3DPD=A3'*Dinv'*P*Dinv; + F1=-(Q+A3DPD*A3)\(A3DPD*A1); %component of (26) + H1=Dinv*(A1+A3*F1); %component of (27) [rcode,NQ]=CheckConvergence([H1;F1]-[H10;F10],iter,solve_maxit,discretion_tol); if rcode @@ -97,16 +137,17 @@ while 1 F10=F1; end +%check if successful retcode = 0; switch rcode case 3 % nan retcode=63; retcode(2)=10000; if verbose - disp([mfilename,':: NAN elements in the solution']) + disp([mfilename,':: NaN elements in the solution']) end case 2% maxiter - retcode = 61 + retcode = 61; if verbose disp([mfilename,':: Maximum Number of Iterations reached']) end @@ -125,8 +166,8 @@ if retcode(1) H=[]; G=[]; else - F2=-(Q+A3DPD*A3)\(A3DPD*A5); - H2=Dinv*(A5+A3*F2); + F2=-(Q+A3DPD*A3)\(A3DPD*A5); %equation (29) + H2=Dinv*(A5+A3*F2); %equation (31) H=zeros(endo_nbr+AuxiliaryVariables_nbr); G=zeros(endo_nbr+AuxiliaryVariables_nbr,exo_nbr); H(endo_augm_id,endo_augm_id)=H1; @@ -159,6 +200,7 @@ end end function [A00,A11,A22,A33,A44,A55,WW,Q,endo_nbr,exo_nbr,aux,endo_augm_id]=GetDennisMatrices(AAlag,AA0,AAlead,BB,bigw,instr_id) +%get the matrices to use the Dennis (2007) algorithm [eq_nbr,endo_nbr]=size(AAlag); exo_nbr=size(BB,2); y=setdiff(1:endo_nbr,instr_id); @@ -211,7 +253,7 @@ end function v = SylvesterHessenbergSchur(d,g,h) % -% DSYLHS Solves a discrete time sylvester equation using the +% DSYLHS Solves a discrete time sylvester equation using the % Hessenberg-Schur algorithm % % v = DSYLHS(g,d,h) computes the matrix v that satisfies the From 9d90a204f76de451b6129b01496a377d837b91e4 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Sun, 30 Aug 2015 11:16:26 +0200 Subject: [PATCH 3/7] Clean up evaluate_planner_objective.m --- matlab/evaluate_planner_objective.m | 19 ++++--------------- 1 file changed, 4 insertions(+), 15 deletions(-) diff --git a/matlab/evaluate_planner_objective.m b/matlab/evaluate_planner_objective.m index e6f80f9bd..69f24c6e5 100644 --- a/matlab/evaluate_planner_objective.m +++ b/matlab/evaluate_planner_objective.m @@ -11,7 +11,7 @@ function planner_objective_value = evaluate_planner_objective(M,options,oo) % SPECIAL REQUIREMENTS % none -% Copyright (C) 2007-2012 Dynare Team +% Copyright (C) 2007-2015 Dynare Team % % This file is part of Dynare. % @@ -29,19 +29,10 @@ function planner_objective_value = evaluate_planner_objective(M,options,oo) % along with Dynare. If not, see . dr = oo.dr; -endo_nbr = M.endo_nbr; exo_nbr = M.exo_nbr; nstatic = M.nstatic; nspred = M.nspred; -lead_lag_incidence = M.lead_lag_incidence; beta = get_optimal_policy_discount_factor(M.params,M.param_names); -if options.ramsey_policy - i_org = (1:M.orig_endo_nbr)'; -else - i_org = (1:M.endo_nbr)'; -end -ipred = find(lead_lag_incidence(M.maximum_lag,:))'; -order_var = dr.order_var; Gy = dr.ghx(nstatic+(1:nspred),:); Gu = dr.ghu(nstatic+(1:nspred),:); @@ -51,11 +42,9 @@ gu(dr.order_var,:) = dr.ghu; ys = oo.dr.ys; -u = oo.exo_simul(1,:)'; - [U,Uy,Uyy] = feval([M.fname '_objective_static'],ys,zeros(1,exo_nbr), ... M.params); - +%second order terms Uyy = full(Uyy); [Uyygygy, err] = A_times_B_kronecker_C(Uyy,gy,gy,options.threads.kronecker.A_times_B_kronecker_C); @@ -65,7 +54,7 @@ mexErrCheck('A_times_B_kronecker_C', err); [Uyygygu, err] = 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); +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)); @@ -75,7 +64,7 @@ mexErrCheck('A_times_B_kronecker_C', err); mexErrCheck('A_times_B_kronecker_C', err); Wuu = Uyygugu+beta*Wyygugu; Wyu = Uyygygu+beta*Wyygygu; -Wss = beta*Wuu*M.Sigma_e(:)/(1-beta); +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; From 3b64c37cb30ec2e3171205dad34a16daa08432f6 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Sun, 30 Aug 2015 11:38:16 +0200 Subject: [PATCH 4/7] Update information in M_ if discretionary_policy_1.m changes lead_lag_incidence --- matlab/discretionary_policy_1.m | 38 +++++++++++++++++++++++++++------ 1 file changed, 32 insertions(+), 6 deletions(-) diff --git a/matlab/discretionary_policy_1.m b/matlab/discretionary_policy_1.m index 4497bd0b3..236e6ec0c 100644 --- a/matlab/discretionary_policy_1.m +++ b/matlab/discretionary_policy_1.m @@ -1,6 +1,6 @@ function [dr,ys,info]=discretionary_policy_1(oo_,Instruments) -% Copyright (C) 2007-2012 Dynare Team +% Copyright (C) 2007-2015 Dynare Team % % This file is part of Dynare. % @@ -88,6 +88,7 @@ end eq_nbr= size(jacobia_,1); instr_nbr=endo_nbr-eq_nbr; + instr_id=nan(instr_nbr,1); for j=1:instr_nbr vj=deblank(Instruments(j,:)); @@ -126,24 +127,49 @@ if info dr=[]; return else - Hold=H; + Hold=H; %save previous solution % Hold=[]; use this line if persistent command is not used. end % update the following elements LLI=lead_lag_incidence; -LLI(MaxLag,:)=any(H); +LLI(MaxLag,:)=any(H); %check if variable drops out in solution LLI=LLI'; tmp=find(LLI); -LLI(tmp)=1:numel(tmp); +LLI(tmp)=1:numel(tmp); %renumber -M_.lead_lag_incidence = LLI'; +M_.lead_lag_incidence = LLI'; %update lead_lag_incidence + +%update info in M_ +max_lag = M_.maximum_endo_lag; +endo_nbr = M_.endo_nbr; +lead_lag_incidence = M_.lead_lag_incidence; + +fwrd_var = find(lead_lag_incidence(max_lag+2:end,:))'; +if max_lag > 0 + pred_var = find(lead_lag_incidence(1,:))'; + both_var = intersect(pred_var,fwrd_var); + pred_var = setdiff(pred_var,both_var); + fwrd_var = setdiff(fwrd_var,both_var); + stat_var = setdiff([1:endo_nbr]',union(union(pred_var,both_var),fwrd_var)); % static variables +else + pred_var = []; + both_var = []; + stat_var = setdiff([1:endo_nbr]',fwrd_var); +end +M_.nstatic=length(stat_var); +M_.nfwrd=length(fwrd_var); +M_.npred=length(pred_var); +M_.nboth=length(both_var); +M_.nspred=M_.npred+M_.nboth; +M_.nsfwrd=M_.nfwrd+M_.nboth; +M_.ndynamic=M_.endo_nbr-M_.nstatic; % set the state dr=oo_.dr; dr.ys =zeros(endo_nbr,1); -dr=set_state_space(dr,M_,options_); +dr=set_state_space(dr,M_,options_); %relies on M_.lead_lag_incidence being updated order_var=dr.order_var; T=H(order_var,order_var); From 9e504a12938b852611f2843024d9f2a1f7b0d0e7 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Sun, 30 Aug 2015 11:43:54 +0200 Subject: [PATCH 5/7] discretionary_policy: Add check whether equation number is consistent with declared number of instruments Also fixes unit test that violated this. --- matlab/discretionary_policy_1.m | 8 +++++++- tests/discretionary_policy/dennis_1.mod | 5 ++--- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/matlab/discretionary_policy_1.m b/matlab/discretionary_policy_1.m index 236e6ec0c..038554618 100644 --- a/matlab/discretionary_policy_1.m +++ b/matlab/discretionary_policy_1.m @@ -88,7 +88,13 @@ end eq_nbr= size(jacobia_,1); instr_nbr=endo_nbr-eq_nbr; - +if instr_nbr==0 + error('discretionary_policy:: There are no available instruments, because the model has as many equations as variables.') +end +if size(Instruments,1)~= instr_nbr + error('discretionary_policy:: There are more declared instruments than omitted equations.') +end + instr_id=nan(instr_nbr,1); for j=1:instr_nbr vj=deblank(Instruments(j,:)); diff --git a/tests/discretionary_policy/dennis_1.mod b/tests/discretionary_policy/dennis_1.mod index f212c9622..0dd883dc3 100644 --- a/tests/discretionary_policy/dennis_1.mod +++ b/tests/discretionary_policy/dennis_1.mod @@ -19,7 +19,7 @@ y = y(+1) -(omega/sigma)*(i-pi(+1))+g; pi = beta*pi(+1)+kappa*y+u; pi_c = pi+(alpha/(1-alpha))*(q-q(-1)); q = q(+1)-(1-alpha)*(i-pi(+1))+(1-alpha)*e; -i = 1.5*pi; +% i = 1.5*pi; end; shocks; @@ -29,6 +29,5 @@ var e; stderr 1; end; planner_objective pi_c^2 + y^2; - -discretionary_policy(instruments=(i),irf=0,qz_criterium=0.999999); +discretionary_policy(instruments=(i),irf=0,planner_discount=beta); From b592987fb4814ab085a059bc022e2f9149a55572 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Sun, 30 Aug 2015 11:56:26 +0200 Subject: [PATCH 6/7] Base evaluation of linearity at steady state for exogenous variables --- matlab/discretionary_policy_1.m | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/matlab/discretionary_policy_1.m b/matlab/discretionary_policy_1.m index 038554618..fd0416304 100644 --- a/matlab/discretionary_policy_1.m +++ b/matlab/discretionary_policy_1.m @@ -78,7 +78,8 @@ it_ = MaxLag + 1 ; if exo_nbr == 0 oo_.exo_steady_state = [] ; end -[junk,jacobia_] = feval([M_.fname '_dynamic'],z, [oo_.exo_simul ... + +[junk,jacobia_] = feval([M_.fname '_dynamic'],z, [zeros(size(oo_.exo_simul)) ... oo_.exo_det_simul], M_.params, zeros(endo_nbr,1), it_); if any(junk~=0) error(['discretionary_policy: the model must be written in deviation ' ... From 752bc543b4a312ee64c97487cea4386a3a9306e7 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Sun, 30 Aug 2015 11:58:37 +0200 Subject: [PATCH 7/7] Add unit test for correctness of discretionary policy --- tests/Makefile.am | 1 + .../discretionary_policy/Gali_discretion.mod | 149 ++++++++++++++++++ 2 files changed, 150 insertions(+) create mode 100644 tests/discretionary_policy/Gali_discretion.mod diff --git a/tests/Makefile.am b/tests/Makefile.am index a639bf838..9afa4a414 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -43,6 +43,7 @@ MODFILES = \ optimal_policy/Ramsey/ramsey_ex_aux.mod \ optimal_policy/RamseyConstraints/test1.mod \ discretionary_policy/dennis_1.mod \ + discretionary_policy/Gali_discretion.mod \ initval_file/ramst_initval_file.mod \ ramst_normcdf_and_friends.mod \ ramst_vec.mod \ diff --git a/tests/discretionary_policy/Gali_discretion.mod b/tests/discretionary_policy/Gali_discretion.mod new file mode 100644 index 000000000..4324a0f47 --- /dev/null +++ b/tests/discretionary_policy/Gali_discretion.mod @@ -0,0 +1,149 @@ +/* + * This file implements the baseline New Keynesian model of Jordi Galí (2008): Monetary Policy, Inflation, + * and the Business Cycle, Princeton University Press, Chapter 5 + * + * This implementation was written by Johannes Pfeifer. + * + * Please note that the following copyright notice only applies to this Dynare + * implementation of the model. + */ + +/* + * Copyright (C) 2013-15 Johannes Pfeifer + * + * This 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. + * + * It 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. + * + * For a copy of the GNU General Public License, + * see . + */ + + +var pi + y_gap + r_e + y_e + r_nat + i + u + a + p + ; + +varexo eps_a + eps_u; + +parameters alppha + betta + rho_a + rho_u + siggma + phi + phi_y + eta + epsilon + theta + ; +%---------------------------------------------------------------- +% Parametrization, p. 52 +%---------------------------------------------------------------- +siggma = 1; +phi=1; +phi_y = .5/4; +theta=2/3; +rho_u = 0; +rho_a = 0.9; +betta = 0.99; +eta =4; +alppha=1/3; +epsilon=6; + + + +%---------------------------------------------------------------- +% First Order Conditions +%---------------------------------------------------------------- + +model(linear); +//Composite parameters +#Omega=(1-alppha)/(1-alppha+alppha*epsilon); //defined on page 47 +#psi_n_ya=(1+phi)/(siggma*(1-alppha)+phi+alppha); //defined on page 48 +#lambda=(1-theta)*(1-betta*theta)/theta*Omega; //defined on page 47 +#kappa=lambda*(siggma+(phi+alppha)/(1-alppha)); //defined on page 49 +#alpha_x=kappa/epsilon; //defined on page 96 +#phi_pi=(1-rho_u)*kappa*siggma/(alpha_x)+rho_u; //defined on page 101 + +r_e=siggma*(y_e(+1)-y_e); +y_e=psi_n_ya*a; +pi=betta*pi(+1)+kappa*y_gap + u; +y_gap=-1/siggma*(i-pi(+1)-r_e)+y_gap(+1); +//3. Interest Rate Rule eq. (25) +% i=r_e+phi_pi*pi; + +r_nat=siggma*psi_n_ya*(a(+1)-a); +u=rho_u*u(-1)+eps_u; +a=rho_a*a(-1)+eps_a; + +pi=p-p(-1); +end; + +%---------------------------------------------------------------- +% define shock variances +%--------------------------------------------------------------- + + +shocks; +var eps_u = 1; +end; + +planner_objective pi^2 +(((1-theta)*(1-betta*theta)/theta*((1-alppha)/(1-alppha+alppha*epsilon)))*(siggma+(phi+alppha)/(1-alppha)))/epsilon*y_gap^2; + +discretionary_policy(instruments=(i),irf=20,planner_discount=betta,discretionary_tol=1e-12) y_gap pi p u; + +verbatim; +%% Check correctness +Omega=(1-alppha)/(1-alppha+alppha*epsilon); %defined on page 47 +lambda=(1-theta)*(1-betta*theta)/theta*Omega; %defined on page 47 +kappa=lambda*(siggma+(phi+alppha)/(1-alppha)); %defined on page 49 +alpha_x=kappa/epsilon; %defined on page 96 +Psi=1/(kappa^2+alpha_x*(1-betta*rho_u)); %defined on page 99 +Psi_i=Psi*(kappa*siggma*(1-rho_u)+alpha_x*rho_u); %defined on page 101 +phi_pi=(1-rho_u)*kappa*siggma/(alpha_x)+rho_u; %defined on page 101 + +%Compute theoretical solution +var_pi_theoretical=(alpha_x*Psi)^2; %equation (6), p.99 +var_y_gap_theoretical=(-kappa*Psi)^2; %equation (7), p.99 + +pi_pos=strmatch('pi',var_list_,'exact'); +y_gap_pos=strmatch('y_gap',var_list_,'exact'); +if abs(oo_.var(pi_pos,pi_pos)-var_pi_theoretical)>1e-10 || abs(oo_.var(y_gap_pos,y_gap_pos)-var_y_gap_theoretical)>1e-10 + error('Variances under optimal policy are wrong') +end + +%Compute theoretical objective function +V=betta/(1-betta)*(var_pi_theoretical+alpha_x*var_y_gap_theoretical); %evaluate at steady state in first period + +if abs(V-oo_.planner_objective_value)>1e-10 + error('Computed welfare deviates from theoretical welfare') +end +end; + +%% repeat exercise with initial shock of 1 to check whether planner objective is correctly specified +initval; + eps_u = 1; +end; + +%Compute theoretical objective function +V=var_pi_theoretical+alpha_x*var_y_gap_theoretical+ betta/(1-betta)*(var_pi_theoretical+alpha_x*var_y_gap_theoretical); %evaluate at steady state in first period + +discretionary_policy(instruments=(i),irf=20,discretionary_tol=1e-12,planner_discount=betta) y_gap pi p u; +if abs(V-oo_.planner_objective_value)>1e-10 + error('Computed welfare deviates from theoretical welfare') +end