discretionary_policy_engine

 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.
 We use the Dennis 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.

0001 function [H,G,retcode]=discretionary_policy_engine(AAlag,AA0,AAlead,BB,bigw,instr_id,beta,solve_maxit,discretion_tol,qz_criterium,H00,verbose)
0002 
0003 % Solves the discretionary problem for a model of the form:
0004 % AAlag*yy_{t-1}+AA0*yy_t+AAlead*yy_{t+1}+BB*e=0, with W the weight on the
0005 % variables in vector y_t and instr_id is the location of the instruments
0006 % in the yy_t vector.
0007 % We use the Dennis algorithm and so we need to re-write the model in the
0008 % 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
0009 % weight on the y_t vector and Q the weight on the x_t vector of
0010 % instruments.
0011 
0012 % Copyright (C) 2007-2011 Dynare Team
0013 %
0014 % This file is part of Dynare.
0015 %
0016 % Dynare is free software: you can redistribute it and/or modify
0017 % it under the terms of the GNU General Public License as published by
0018 % the Free Software Foundation, either version 3 of the License, or
0019 % (at your option) any later version.
0020 %
0021 % Dynare is distributed in the hope that it will be useful,
0022 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0023 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0024 % GNU General Public License for more details.
0025 %
0026 % You should have received a copy of the GNU General Public License
0027 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0028 
0029 if nargin<12
0030     verbose=0;
0031     if nargin<11
0032         H00=[];
0033         if nargin<10
0034             qz_criterium=1.000001;
0035             if nargin<9
0036                 discretion_tol=sqrt(eps);
0037                 if nargin<8
0038                     solve_maxit=3000;
0039                     if nargin<7
0040                         beta=.99;
0041                         if nargin<6
0042                             error([mfilename,':: Insufficient number of input arguments'])
0043                         elseif nargin>12
0044                             error([mfilename,':: Number of input arguments cannot exceed 12'])
0045                         end
0046                     end
0047                 end
0048             end
0049         end
0050     end
0051 end
0052 
0053 [A0,A1,A2,A3,A4,A5,W,Q,endo_nbr,exo_nbr,aux,endo_augm_id]=GetDennisMatrices(AAlag,AA0,AAlead,BB,bigw,instr_id);
0054 % aux is a logical index of the instruments which appear with lags in the
0055 % model. Their location in the state vector is instr_id(aux)
0056 % endo_augm_id is index (not logical) of locations of the augmented vector
0057 % of non-instrumental variables
0058 
0059 AuxiliaryVariables_nbr=sum(aux);
0060 H0=zeros(endo_nbr+AuxiliaryVariables_nbr);
0061 if ~isempty(H00)
0062     H0(1:endo_nbr,1:endo_nbr)=H00;clear H00
0063 end
0064 
0065 H10=H0(endo_augm_id,endo_augm_id);
0066 F10=H0(instr_id,endo_augm_id);
0067 
0068 iter=0;
0069 H1=H10;
0070 F1=F10;
0071 while 1
0072     iter=iter+1;
0073     P=SylvesterDoubling(W+beta*F1'*Q*F1,beta*H1',H1,discretion_tol,solve_maxit);
0074     if any(any(isnan(P)))
0075         P=SylvesterHessenbergSchur(W+beta*F1'*Q*F1,beta*H1',H1);
0076         if any(any(isnan(P)))
0077             retcode=2;
0078             return
0079         end
0080     end
0081     D=A0-A2*H1-A4*F1;
0082     Dinv=inv(D);
0083     A3DPD=A3'*Dinv'*P*Dinv;
0084     F1=-(Q+A3DPD*A3)\(A3DPD*A1);
0085     H1=Dinv*(A1+A3*F1);
0086     
0087     [rcode,NQ]=CheckConvergence([H1;F1]-[H10;F10],iter,solve_maxit,discretion_tol);
0088     if rcode
0089         break
0090     else
0091         if verbose
0092             disp(NQ)
0093         end
0094     end
0095     H10=H1;
0096     F10=F1;
0097 end
0098 
0099 retcode = 0;
0100 switch rcode
0101   case 3 % nan
0102     retcode=63;
0103     retcode(2)=10000;
0104     if verbose
0105         disp([mfilename,':: NAN elements in the solution'])
0106     end
0107   case 2% maxiter
0108     retcode = 61
0109     if verbose
0110         disp([mfilename,':: Maximum Number of Iterations reached'])
0111     end
0112   case 1
0113     BadEig=max(abs(eig(H1)))>qz_criterium;
0114     if BadEig
0115         retcode=62;
0116         retcode(2)=100*max(abs(eig(H1)));
0117         if verbose
0118             disp([mfilename,':: Some eigenvalues greater than qz_criterium, Model potentially unstable'])
0119         end
0120     end
0121 end
0122 
0123 if retcode(1)
0124     H=[];
0125     G=[];
0126 else
0127     F2=-(Q+A3DPD*A3)\(A3DPD*A5);
0128     H2=Dinv*(A5+A3*F2);
0129     H=zeros(endo_nbr+AuxiliaryVariables_nbr);
0130     G=zeros(endo_nbr+AuxiliaryVariables_nbr,exo_nbr);
0131     H(endo_augm_id,endo_augm_id)=H1;
0132     H(instr_id,endo_augm_id)=F1;
0133     G(endo_augm_id,:)=H2;
0134     G(instr_id,:)=F2;
0135     
0136     % Account for auxilliary variables
0137     H(:,instr_id(aux))=H(:,end-(AuxiliaryVariables_nbr-1:-1:0));
0138     H=H(1:endo_nbr,1:endo_nbr);
0139     G=G(1:endo_nbr,:);
0140 end
0141 
0142 end
0143 
0144 
0145 function [rcode,NQ]=CheckConvergence(Q,iter,MaxIter,crit)
0146 
0147 NQ=max(max(abs(Q)));% norm(Q); seems too costly
0148 if isnan(NQ)
0149     rcode=3;
0150 elseif iter>MaxIter;
0151     rcode=2;
0152 elseif NQ<crit
0153     rcode=1;
0154 else
0155     rcode=0;
0156 end
0157 
0158 end
0159 
0160 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)
0161 [eq_nbr,endo_nbr]=size(AAlag);
0162 exo_nbr=size(BB,2);
0163 y=setdiff(1:endo_nbr,instr_id);
0164 instr_nbr=numel(instr_id);
0165 
0166 A0=AA0(:,y);
0167 A1=-AAlag(:,y);
0168 A2=-AAlead(:,y);
0169 A3=-AA0(:,instr_id);
0170 A4=-AAlead(:,instr_id);
0171 A5=-BB;
0172 W=bigw(y,y);
0173 Q=bigw(instr_id,instr_id);
0174 % Adjust for possible lags in instruments by creating auxiliary equations
0175 A6=-AAlag(:,instr_id);
0176 aux=any(A6);
0177 AuxiliaryVariables_nbr=sum(aux);
0178 ny=eq_nbr;
0179 m=eq_nbr+AuxiliaryVariables_nbr;
0180 A00=zeros(m);A00(1:ny,1:ny)=A0;A00(ny+1:end,ny+1:end)=eye(AuxiliaryVariables_nbr);
0181 A11=zeros(m);A11(1:ny,1:ny)=A1;A11(1:ny,ny+1:end)=A6(:,aux);
0182 A22=zeros(m);A22(1:ny,1:ny)=A2;
0183 A33=zeros(m,instr_nbr);A33(1:ny,1:end)=A3;A33(ny+1:end,aux)=eye(AuxiliaryVariables_nbr);
0184 A44=zeros(m,instr_nbr);A44(1:ny,1:end)=A4;
0185 A55=zeros(m,exo_nbr);A55(1:ny,1:end)=A5;
0186 WW=zeros(m);WW(1:ny,1:ny)=W;
0187 endo_augm_id=setdiff(1:endo_nbr+AuxiliaryVariables_nbr,instr_id);
0188 
0189 end
0190 
0191 function v= SylvesterDoubling (d,g,h,tol,maxit)
0192 
0193 % DOUBLES  Solves a Sylvester equation using doubling
0194 %
0195 %  [v,info] = doubles (g,d,h,tol,maxit)  uses a doubling algorithm
0196 %   to solve the  Sylvester equation v  = d + g v h
0197 
0198 v = d;
0199 for i =1:maxit,
0200     vadd = g*v*h;
0201     v = v+vadd;
0202     if norm (vadd,1) <= (tol*norm(v,1))
0203         break;
0204     end
0205     g = g*g;
0206     h = h*h;
0207 end
0208 
0209 end
0210 
0211 function v = SylvesterHessenbergSchur(d,g,h)
0212 %
0213 % DSYLHS  Solves a discrete time sylvester equation     using the
0214 % Hessenberg-Schur algorithm
0215 %
0216 % v = DSYLHS(g,d,h) computes the matrix v that satisfies the
0217 % discrete time sylvester equation
0218 %
0219 %           v = d + g'vh
0220 
0221 if size(g,1) >= size(h,1)
0222     [u,gbarp] = hess(g');
0223     [t,hbar] = schur(h);
0224     [vbar] = sylvest_private(gbarp,u'*d*t,hbar,1e-15);
0225     v = u*vbar*t';
0226 else
0227     [u,gbar] = schur(g);
0228     [t,hbarp] = hess(h');
0229     [vbar] = sylvest_private(hbarp,t'*d'*u,gbar,1e-15);
0230     v = u*vbar'*t';
0231 end
0232 
0233 end
0234 
0235 
0236 function v = sylvest_private(g,d,h,tol)
0237 %
0238 % SYLVEST  Solves a Sylvester equation
0239 %
0240 %  solves the Sylvester equation
0241 %         v = d + g v h
0242 %  for v where both g and h must be  upper block triangular.
0243 %  The output info is zero on a successful return.
0244 %  The input tol indicates when an element of g or h should be considered
0245 %  zero.
0246 
0247 [m,n] = size(d);
0248 v = zeros(m,n);
0249 w = eye(m);
0250 i = 1;
0251 temp = [];
0252 
0253 %First handle the i = 1 case outside the loop
0254 
0255 if i< n,
0256     if abs(h(i+1,i)) < tol,
0257         v(:,i)= (w - g*h(i,i))\d(:,i);
0258         i = i+1;
0259     else
0260         A = [w-g*h(i,i)     (-g*h(i+1,i));...
0261             -g*h(i,i+1)    w-g*h(i+1,i+1)];
0262         C = [d(:,i); d(:,i+1)];
0263         X = A\C;
0264         v(:,i) = X(1:m,:);
0265         v(:,i+1) = X(m+1:2*m,  :);
0266         i = i+2;
0267     end
0268 end
0269 
0270 %Handle the rest of the matrix with the possible exception of i=n
0271 
0272 while i<n,
0273     b= i-1;
0274     temp = [temp g*v(:,size(temp,2)+1:b)]; %#ok<AGROW>
0275     if abs(h(i+1,i)) < tol,
0276         v(:,i) = (w - g*h(i,i))\(d(:,i) + temp*h(1:b,i));
0277         i = i+1;
0278     else
0279         A = [w - g*h(i,i)    (-g*h(i+1,i));   ...
0280             -g*h(i,i+1)    w - g*h(i+1,i+1)];
0281         C = [d(:,i) + temp*h(1:b,i);         ...
0282             d(:,i+1) + temp*h(1:b,i+1)];
0283         X = A\C;
0284         v(:,i) = X(1:m,:);
0285         v(:,i+1) = X(m+1:2*m, :);
0286         i = i+2;
0287     end
0288 end
0289 
0290 %Handle the i = n case if i=n was not in a 2-2 block
0291 
0292 if i==n,
0293     b = i-1;
0294     temp = [temp g*v(:,size(temp,2)+1:b)];
0295     v(:,i) = (w-g*h(i,i))\(d(:,i) + temp*h(1:b,i));
0296 end
0297 
0298 end