1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273
|
function [G1pi,C,impact,nmat,TT1,TT2,gev,eu, DD, E2, E5, GAMMA, FL_RANK ]=PI_gensys(a0,a1,a2,a3,c,PSI,NX,NETA,FL_RANK,M_,options_)
% System given as
% a0*E_t[y(t+1])+a1*y(t)=a2*y(t-1)+c+psi*eps(t)
% with z an exogenous variable process.
% Returned system is
% [s(t)' x(t)' E_t x(t+1)']'=G1pi [s(t-1)' x(t-1)' x(t)]'+C+impact*eps(t),
% and (a) the matrix nmat satisfying nmat*E_t z(t)+ E_t x(t+1)=0
% (b) matrices TT1, TT2 that relate y(t) to these states: y(t)=[TT1 TT2][s(t)' x(t)']'.
% Note that the dimension of the state vector = dim(a0)+NO_FL_EQS
%
% If div is omitted from argument list, a div>1 is calculated.
% eu(1)=1 for existence, eu(2)=1 for uniqueness. eu(1)=-1 for
% existence only with not-s.c. z; eu=[-2,-2] for coincident zeros.
% Based on
% Christopher A. Sims
% Corrected 10/28/96 by CAS
% Copyright (C) 1996-2009 Christopher Sims
% Copyright (C) 2010-2017 Dynare Team
%
% This file is part of Dynare.
%
% Dynare 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.
%
% Dynare 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.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
lastwarn('','');
global lq_instruments;
eu=[0;0];C=c;
realsmall=1e-6;
fixdiv=(nargin==6);
n=size(a0,1);
DD=[];E2=[]; E5=0; GAMMA=[];
%
% Find SVD of a0, and create partitions of U, S and V
%
[U0,S0,V0] = svd(a0);
FL_RANK=rank(S0);
U01=U0(1:n,1:FL_RANK);
U02=U0(1:n,FL_RANK+1:n);
V01=V0(1:n,1:FL_RANK);
V02=V0(1:n,FL_RANK+1:n);
S01=S0(1:FL_RANK,1:FL_RANK);
C1=U02'*a1*V01;
C2=U02'*a1*V02;
C3=U02'*a2*V01;
C4=U02'*a2*V02;
C5=U02'*PSI;
Sinv=eye(FL_RANK);
for i=1:FL_RANK
Sinv(i,i)=1/S01(i,i);
end
F1=Sinv*U01'*a1*V01;
F2=Sinv*U01'*a1*V02;
F3=Sinv*U01'*a2*V01;
F4=Sinv*U01'*a2*V02;
F5=Sinv*U01'*PSI;
singular=0;
warning('', '');
try
if rcond(C2) < 1e-8
singular=1;
else
warning('off','MATLAB:nearlySingularMatrix');
warning('off','MATLAB:singularMatrix');
UAVinv=inv(C2); % i.e. inv(U02'*a1*V02)
[LastWarningTxt LastWarningID]=lastwarn;
if any(any(isinf(UAVinv)))==1
singular=1;
end
end
if singular == 1 || strcmp('MATLAB:nearlySingularMatrix',LastWarningID) == 1 || ...
strcmp('MATLAB:illConditionedMatrix',LastWarningID)==1 || ...
strcmp('MATLAB:singularMatrix',LastWarningID)==1
[C1,C2,C3,C4, C5, F1, F2, F3, F4, F5, M1, M2, UAVinv, FL_RANK, V01, V02] = PI_gensys_singularC(C1,C2,C3,C4, C5, F1, F2, F3, F4, F5, V01, V02, 0);
end
warning('on','MATLAB:singularMatrix');
warning('on','MATLAB:nearlySingularMatrix');
if (any(any(isinf(UAVinv))) || any(any(isnan(UAVinv))))
if(options_.ACES_solver==1)
disp('ERROR! saving PI_gensys_data_dump');
save PI_gensys_data_dump
error('PI_gensys: Inversion of poss. zero matrix UAVinv=inv(U02''*a1*V02)!');
else
warning('PI_gensys: Evading inversion of zero matrix UAVinv=inv(U02''*a1*V02)!');
eu=[0,0];
return
end
end
catch
errmsg=lasterror;
warning(['error callig PI_gensys_singularC: ' errmsg.message ],'errmsg.identifier');
%error('errcode',['error callig PI_gensys_singularC: ' errmsg.message ]);
end
%
% Define TT1, TT2
%
TT1=V02;
TT2=V01;
%
%Set up the system matrix for the variables s(t)=V02*Y(t), x(t)=V01*Y(t) and E_t x(t+1)
% and define z(t)'=[s(t)' x(t)]
%
%UAVinv=inv(U02'*a1*V02);
FF=F2; %=Sinv*U01'*a1*V02;
G11=UAVinv*C4;
G12=UAVinv*C3;
G13=-UAVinv*C1;
G31=-FF*G11+F4; % +Sinv*U01'*a2*V02;
G32=-FF*G12+F3; % +Sinv*U01'*a2*V01;
G33=-FF*G13-F1; % -Sinv*U01'*a1*V01;
P1=UAVinv*C5; % *U02'*PSI;
P3=-FF*P1+F5; % + Sinv*U01'*PSI; % This should equal 0
G21=zeros(FL_RANK,(n-FL_RANK));
G22=zeros(FL_RANK,FL_RANK);
G23=eye(FL_RANK);
%H2=zeros(FL_RANK,NX);
num_inst=0;
% New Definitions
Ze11=zeros(NX,NX);
Ze12=zeros(NX,(n-FL_RANK));
Ze134=zeros(NX,FL_RANK);
Ze31=zeros(FL_RANK,NX);
% End of New Definitions
%
% System matrix is called 'G1pi'; Shock matrix is called 'impact'
%
G1pi=[Ze11 Ze12 Ze134 Ze134; P1 G11 G12 G13; Ze31 G21 G22 G23; P3 G31 G32 G33];
impact=[eye(NX,NX); zeros(n+FL_RANK,NX)];
if(options_.ACES_solver==1)
if isfield(lq_instruments,'names')
num_inst=size(lq_instruments.names,1);
if num_inst>0
i_var=lq_instruments.inst_var_indices;
N1=UAVinv*U02'*lq_instruments.B1;
N3=-FF*N1+Sinv*U01'*lq_instruments.B1;
else
error('WARNING: There are no instrumnets for ACES!');
end
lq_instruments.N1=N1;
lq_instruments.N3=N3;
else
error('WARNING: There are no instrumnets for ACES!');
end
E3=V02*[P1 G11 G12 G13];
E3= E3+ [zeros(size(V01,1),size(E3,2)-size(V01,2)) V01];
E2=-E3;
E5=-V02*N1;
DD=[zeros(NX,size(N1,2));N1; zeros(FL_RANK,size(N1,2));N3];
II=eye(num_inst);
GAMMA=[ E3 -E5 %zeros(size(E3,1),num_inst);
zeros(num_inst,size(E3,2)), II;
];
eu =[1; 1], nmat=[], gev=[];
return % do not check B&K compliancy
end
G0pi=eye(n+FL_RANK+NX);
try
if isoctave
% Need to force QZ complex on Octave (otherwise it returns the real one)
[a b q z]=qz(complex(G0pi),complex(G1pi));
else
[a b q z]=qz(G0pi,G1pi);
end
catch
try
lerror=lasterror;
disp(['PI_Gensys: ' lerror.message]);
if 0==strcmp('MATLAB:qz:matrixWithNaNInf',lerror.identifier)
disp '** Unexpected Error PI_Gensys:qz: ** :';
button=questdlg('Continue Y/N?','Unexpected Error in qz','No','Yes','Yes');
switch button
case 'No'
error ('Terminated')
%case 'Yes'
end
end
G1pi=[];impact=[];nmat=[]; gev=[];
eu=[-2;-2];
return
catch
disp '** Unexpected Error in qz ** :';
disp lerror.message;
button=questdlg('Continue Y/N?','Unexpected Error in qz','No','Yes','Yes');
switch button
case 'No'
error ('Terminated')
case 'Yes'
G1pi=[];impact=[];nmat=[]; gev=[];
eu=[-2;-2];
return
end
end
end
if ~fixdiv, div=1.01; end
nunstab=0;
zxz=0;
nn=size(a,1);
for i=1:nn
% ------------------div calc------------
if ~fixdiv
if abs(a(i,i)) > 0
divhat=abs(b(i,i))/abs(a(i,i));
% bug detected by Vasco Curdia and Daria Finocchiaro, 2/25/2004 A root of
% exactly 1.01 and no root between 1 and 1.02, led to div being stuck at 1.01
% and the 1.01 root being misclassified as stable. Changing < to <= below fixes this.
if 1+realsmall<divhat && divhat<=div
div=.5*(1+divhat);
end
end
end
% ----------------------------------------
nunstab=nunstab+(abs(b(i,i))>div*abs(a(i,i)));
if abs(a(i,i))<realsmall && abs(b(i,i))<realsmall
zxz=1;
end
end
div ;
if ~zxz
[a b q z]=qzdiv(div,a,b,q,z);
end
gev=[diag(a) diag(b)];
if zxz
disp('Coincident zeros. Indeterminacy and/or nonexistence.')
eu=[-2;-2];
% correction added 7/29/2003. Otherwise the failure to set output
% arguments leads to an error message and no output (including eu).
nmat=[]; %;gev=[]
return
end
if (FL_RANK ~= nunstab && options_.ACES_solver~=1)
disp(['Number of unstable variables ' num2str(nunstab)]);
disp( ['does not match number of expectational equations ' num2str(FL_RANK)]);
nmat=[];% gev=[];
eu=[-2;-2];
return
end
% New Definitions
z1=z(:,1:n+NX)';
z2=z(:,n+NX+1:n+NX+FL_RANK)';
% New N Matrix by J Pearlman
z12=z2(:,1:n+NX);
z22=z2(:,n+NX+1:n+NX+FL_RANK);
% End of New Definitions
% modified by GP:
nmat=real(inv(z22)*z12);
eu=[1;1];
|