1function [H,G,retcode]=discretionary_policy_engine(AAlag,AA0,AAlead,BB,bigw,instr_id,beta,solve_maxit,discretion_tol,qz_criterium,H00,verbose) 2 3% Solves the discretionary problem for a model of the form: 4% 5% Loss=E_0 sum_{t=0}^{\infty} beta^t [y_t'*W*y+x_t'*Q*x_t] 6% subject to 7% AAlag*yy_{t-1}+AA0*yy_t+AAlead*yy_{t+1}+BB*e=0 8% 9% with W the weight on the variables in vector y_t. 10% 11% The solution takes the form 12% y_t=H*y_{t-1}+G*e_t 13% where H=[H1;F1] and G=[H2;F2]. 14% 15% We use the Dennis (2007, Macroeconomic Dynamics) algorithm and so we need 16% to re-write the model in the form 17% A0*y_t=A1*y_{t-1}+A2*y_{t+1}+A3*x_t+A4*x_{t+1}+A5*e_t, with W the 18% weight on the y_t vector and Q the weight on the x_t vector of 19% instruments. 20% 21% Inputs: 22% AAlag [double] matrix of coefficients on lagged 23% variables 24% AA0 [double] matrix of coefficients on 25% contemporaneous variables 26% AAlead [double] matrix of coefficients on 27% leaded variables 28% BB [double] matrix of coefficients on 29% shocks 30% bigw [double] matrix of coefficients on variables in 31% loss/objective function; stacks [W and Q] 32% instr_id [double] location vector of the instruments in the yy_t vector. 33% beta [scalar] planner discount factor 34% solve_maxit [scalar] maximum number of iterations 35% discretion_tol [scalar] convergence criterion for solution 36% qz_criterium [scalar] tolerance for QZ decomposition 37% H00 38% verbose [scalar] dummy to control verbosity 39% 40% Outputs: 41% H [double] (endo_nbr*endo_nbr) solution matrix for endogenous 42% variables, stacks [H1 and H1] 43% G [double] (endo_nbr*exo_nbr) solution matrix for shocks, stacks [H2 and F2] 44% 45% retcode [scalar] return code 46% 47% Algorithm: 48% Dennis, Richard (2007): Optimal policy in rational expectations models: new solution algorithms, 49% Macroeconomic Dynamics, 11, 3155. 50 51% Copyright (C) 2007-2018 Dynare Team 52% 53% This file is part of Dynare. 54% 55% Dynare is free software: you can redistribute it and/or modify 56% it under the terms of the GNU General Public License as published by 57% the Free Software Foundation, either version 3 of the License, or 58% (at your option) any later version. 59% 60% Dynare is distributed in the hope that it will be useful, 61% but WITHOUT ANY WARRANTY; without even the implied warranty of 62% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 63% GNU General Public License for more details. 64% 65% You should have received a copy of the GNU General Public License 66% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 67 68if nargin<12 69 verbose=0; 70 if nargin<11 71 H00=[]; 72 if nargin<10 73 qz_criterium=1.000001; 74 if nargin<9 75 discretion_tol=sqrt(eps); 76 if nargin<8 77 solve_maxit=3000; 78 if nargin<7 79 beta=.99; 80 if nargin<6 81 error([mfilename,':: Insufficient number of input arguments']) 82 elseif nargin>12 83 error([mfilename,':: Number of input arguments cannot exceed 12']) 84 end 85 end 86 end 87 end 88 end 89 end 90end 91 92[A0,A1,A2,A3,A4,A5,W,Q,endo_nbr,exo_nbr,aux,endo_augm_id]=GetDennisMatrices(AAlag,AA0,AAlead,BB,bigw,instr_id); 93% aux is a logical index of the instruments which appear with lags in the 94% model. Their location in the state vector is instr_id(aux); 95% endo_augm_id is index (not logical) of locations of the augmented vector 96% of non-instrumental variables 97 98AuxiliaryVariables_nbr=sum(aux); 99H0=zeros(endo_nbr+AuxiliaryVariables_nbr); 100if ~isempty(H00) 101 H0(1:endo_nbr,1:endo_nbr)=H00; 102 clear H00 103end 104 105H10=H0(endo_augm_id,endo_augm_id); 106F10=H0(instr_id,endo_augm_id); 107 108iter=0; 109H1=H10; 110F1=F10; 111%solve equations (20) and (22) via fixed point iteration 112while 1 113 iter=iter+1; 114 P=SylvesterDoubling(W+beta*F1'*Q*F1,beta*H1',H1,discretion_tol,solve_maxit); 115 if any(any(isnan(P))) || any(any(isinf(P))) 116 P=SylvesterHessenbergSchur(W+beta*F1'*Q*F1,beta*H1',H1); 117 if any(any(isnan(P))) || any(any(isinf(P))) 118 retcode=2; 119 return 120 end 121 end 122 D=A0-A2*H1-A4*F1; %equation (20) 123 Dinv=inv(D); 124 A3DPD=A3'*Dinv'*P*Dinv; 125 F1=-(Q+A3DPD*A3)\(A3DPD*A1); %component of (26) 126 H1=Dinv*(A1+A3*F1); %component of (27) 127 128 [rcode,NQ]=CheckConvergence([H1;F1]-[H10;F10],iter,solve_maxit,discretion_tol); 129 if rcode 130 break 131 else 132 if verbose 133 disp(NQ) 134 end 135 end 136 H10=H1; 137 F10=F1; 138end 139 140%check if successful 141retcode = 0; 142switch rcode 143 case 3 % nan 144 retcode=63; 145 retcode(2)=10000; 146 if verbose 147 disp([mfilename,':: NaN elements in the solution']) 148 end 149 case 2% maxiter 150 retcode = 61; 151 if verbose 152 disp([mfilename,':: Maximum Number of Iterations reached']) 153 end 154 case 1 155 BadEig=max(abs(eig(H1)))>qz_criterium; 156 if BadEig 157 retcode=62; 158 retcode(2)=100*max(abs(eig(H1))); 159 if verbose 160 disp([mfilename,':: Some eigenvalues greater than qz_criterium, Model potentially unstable']) 161 end 162 end 163end 164 165if retcode(1) 166 H=[]; 167 G=[]; 168else 169 F2=-(Q+A3DPD*A3)\(A3DPD*A5); %equation (29) 170 H2=Dinv*(A5+A3*F2); %equation (31) 171 H=zeros(endo_nbr+AuxiliaryVariables_nbr); 172 G=zeros(endo_nbr+AuxiliaryVariables_nbr,exo_nbr); 173 H(endo_augm_id,endo_augm_id)=H1; 174 H(instr_id,endo_augm_id)=F1; 175 G(endo_augm_id,:)=H2; 176 G(instr_id,:)=F2; 177 178 % Account for auxilliary variables 179 H(:,instr_id(aux))=H(:,end-(AuxiliaryVariables_nbr-1:-1:0)); 180 H=H(1:endo_nbr,1:endo_nbr); 181 G=G(1:endo_nbr,:); 182end 183 184end 185 186 187function [rcode,NQ]=CheckConvergence(Q,iter,MaxIter,crit) 188 189NQ=max(max(abs(Q)));% norm(Q); seems too costly 190if isnan(NQ) 191 rcode=3; 192elseif iter>MaxIter; 193 rcode=2; 194elseif NQ<crit 195 rcode=1; 196else 197 rcode=0; 198end 199 200end 201 202function [A00,A11,A22,A33,A44,A55,WW,Q,endo_nbr,exo_nbr,aux,endo_augm_id]=GetDennisMatrices(AAlag,AA0,AAlead,BB,bigw,instr_id) 203%get the matrices to use the Dennis (2007) algorithm 204[eq_nbr,endo_nbr]=size(AAlag); 205exo_nbr=size(BB,2); 206y=setdiff(1:endo_nbr,instr_id); 207instr_nbr=numel(instr_id); 208 209A0=AA0(:,y); 210A1=-AAlag(:,y); 211A2=-AAlead(:,y); 212A3=-AA0(:,instr_id); 213A4=-AAlead(:,instr_id); 214A5=-BB; 215W=bigw(y,y); 216Q=bigw(instr_id,instr_id); 217% Adjust for possible lags in instruments by creating auxiliary equations 218A6=-AAlag(:,instr_id); 219aux=any(A6); 220AuxiliaryVariables_nbr=sum(aux); 221ny=eq_nbr; 222m=eq_nbr+AuxiliaryVariables_nbr; 223A00=zeros(m);A00(1:ny,1:ny)=A0;A00(ny+1:end,ny+1:end)=eye(AuxiliaryVariables_nbr); 224A11=zeros(m);A11(1:ny,1:ny)=A1;A11(1:ny,ny+1:end)=A6(:,aux); 225A22=zeros(m);A22(1:ny,1:ny)=A2; 226A33=zeros(m,instr_nbr);A33(1:ny,1:end)=A3;A33(ny+1:end,aux)=eye(AuxiliaryVariables_nbr); 227A44=zeros(m,instr_nbr);A44(1:ny,1:end)=A4; 228A55=zeros(m,exo_nbr);A55(1:ny,1:end)=A5; 229WW=zeros(m);WW(1:ny,1:ny)=W; 230endo_augm_id=setdiff(1:endo_nbr+AuxiliaryVariables_nbr,instr_id); 231 232end 233 234function v= SylvesterDoubling (d,g,h,tol,maxit) 235 236% DOUBLES Solves a Sylvester equation using doubling 237% 238% [v,info] = doubles (g,d,h,tol,maxit) uses a doubling algorithm 239% to solve the Sylvester equation v = d + g v h 240 241v = d; 242for i =1:maxit 243 vadd = g*v*h; 244 v = v+vadd; 245 if norm (vadd,1) <= (tol*norm(v,1)) 246 break 247 end 248 g = g*g; 249 h = h*h; 250end 251 252end 253 254function v = SylvesterHessenbergSchur(d,g,h) 255% 256% DSYLHS Solves a discrete time sylvester equation using the 257% Hessenberg-Schur algorithm 258% 259% v = DSYLHS(g,d,h) computes the matrix v that satisfies the 260% discrete time sylvester equation 261% 262% v = d + g'vh 263 264if size(g,1) >= size(h,1) 265 [u,gbarp] = hess(g'); 266 [t,hbar] = schur(h); 267 [vbar] = sylvest_private(gbarp,u'*d*t,hbar,1e-15); 268 v = u*vbar*t'; 269else 270 [u,gbar] = schur(g); 271 [t,hbarp] = hess(h'); 272 [vbar] = sylvest_private(hbarp,t'*d'*u,gbar,1e-15); 273 v = u*vbar'*t'; 274end 275 276end 277 278 279function v = sylvest_private(g,d,h,tol) 280% 281% SYLVEST Solves a Sylvester equation 282% 283% solves the Sylvester equation 284% v = d + g v h 285% for v where both g and h must be upper block triangular. 286% The output info is zero on a successful return. 287% The input tol indicates when an element of g or h should be considered 288% zero. 289 290[m,n] = size(d); 291v = zeros(m,n); 292w = eye(m); 293i = 1; 294temp = []; 295 296%First handle the i = 1 case outside the loop 297 298if i< n 299 if abs(h(i+1,i)) < tol 300 v(:,i)= (w - g*h(i,i))\d(:,i); 301 i = i+1; 302 else 303 A = [w-g*h(i,i) (-g*h(i+1,i));... 304 -g*h(i,i+1) w-g*h(i+1,i+1)]; 305 C = [d(:,i); d(:,i+1)]; 306 X = A\C; 307 v(:,i) = X(1:m,:); 308 v(:,i+1) = X(m+1:2*m, :); 309 i = i+2; 310 end 311end 312 313%Handle the rest of the matrix with the possible exception of i=n 314 315while i<n 316 b= i-1; 317 temp = [temp g*v(:,size(temp,2)+1:b)]; %#ok<AGROW> 318 if abs(h(i+1,i)) < tol 319 v(:,i) = (w - g*h(i,i))\(d(:,i) + temp*h(1:b,i)); 320 i = i+1; 321 else 322 A = [w - g*h(i,i) (-g*h(i+1,i)); ... 323 -g*h(i,i+1) w - g*h(i+1,i+1)]; 324 C = [d(:,i) + temp*h(1:b,i); ... 325 d(:,i+1) + temp*h(1:b,i+1)]; 326 X = A\C; 327 v(:,i) = X(1:m,:); 328 v(:,i+1) = X(m+1:2*m, :); 329 i = i+2; 330 end 331end 332 333%Handle the i = n case if i=n was not in a 2-2 block 334 335if i==n 336 b = i-1; 337 temp = [temp g*v(:,size(temp,2)+1:b)]; 338 v(:,i) = (w-g*h(i,i))\(d(:,i) + temp*h(1:b,i)); 339end 340 341end 342