0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020 function [S0, S1, S2, S3, Sy, Sf, Sz, unique] = Smats(C, GAM0, GAM1, PSI, PPI);
0021
0022
0023
0024
0025 eu=[0;0];
0026
0027 realsmall=1e-6;
0028
0029 [n k] = size(PPI);
0030 [n l] = size(PSI);
0031
0032 [A B Q Z VV]=qz(GAM0,GAM1);
0033
0034 div = 1.01;
0035 nunstab=0;
0036 zxz=0;
0037
0038 for i=1:n
0039 if abs(A(i,i)) > 0
0040 divhat=abs(B(i,i))/abs(A(i,i));
0041 if 1+realsmall<divhat && divhat<=div
0042 div=.5*(1+divhat);
0043 end
0044 end
0045 nunstab = nunstab + (abs(B(i,i))>div*abs(A(i,i)));
0046 if abs(A(i,i))<realsmall && abs(B(i,i))<realsmall
0047 zxz=1;
0048 end
0049 end
0050
0051 div ;
0052 nunstab;
0053 m = nunstab;
0054
0055 if ~zxz
0056 [A B Q Z]= qzdiv(div,A,B,Q,Z);
0057 end
0058
0059
0060
0061
0062 gev=[diag(A) diag(B)];
0063
0064 if zxz
0065 disp('Coincident zeros. Indeterminacy and/or nonexistence.')
0066 eu=[-2;-2];
0067
0068
0069
0070
0071
0072
0073 S0 = []; S1=[]; S2=[]; S3=[]; unique=[];
0074 return
0075 end
0076
0077
0078 Q1=Q(1:n-nunstab,:);
0079 Q2=Q(n-nunstab+1:n,:);
0080 Z1 = Z(:,1:n-nunstab)';
0081 Z2 = Z(:,n-nunstab+1:n)';
0082 A22 = A(n-nunstab+1:n,n-nunstab+1:n);
0083 B22 = B(n-nunstab+1:n,n-nunstab+1:n);
0084 etawt=Q2*PPI;
0085
0086
0087
0088 A11 = A(1:n-nunstab,1:n-nunstab);
0089 A12 = A(1:n-nunstab,n-nunstab+1:n);
0090
0091
0092
0093 [U,D,V] = svd(etawt);
0094 md = min(size(D));
0095 bigev = find(diag(D(1:md,1:md))>realsmall);
0096 U1 = U(:,bigev);
0097 V1 = V(:,bigev);
0098 D11 = D(bigev,bigev);
0099 r = size(D11,1);
0100 V2 = V(:,r+1:end);
0101 U2 = U(:,r+1:end);
0102
0103 eu(1) = length(bigev)>=nunstab;
0104
0105 [Uhat,Dhat,Vhat] = svd(Q1*PPI);
0106 md = min(size(Dhat));
0107 bigev = find(diag(Dhat(1:md,1:md))>realsmall);
0108 Uhat1 = Uhat(:,bigev);
0109 Vhat1 = Vhat(:,bigev);
0110 Dhat11 = Dhat(bigev,bigev);
0111 rr = size(Dhat11,1);
0112 Uhat2 = Uhat(:,rr+1:end);
0113 Vhat2 = Vhat(:,rr+1:end);
0114
0115
0116 if eu(1)
0117 disp('Solution Exists')
0118 else
0119 disp('Solution Does Not Exist')
0120 eu
0121 disp('Generalized eigenvalues')
0122
0123
0124
0125 S0 = []; S1=[]; S2=[]; S3=[]; unique=0; Sy=[]; Sf=[]; Sz=[];
0126 return
0127 end
0128
0129 if isempty(V2)
0130 unique=1;
0131 disp('The solution is unique');
0132
0133 else
0134 unique=0;
0135 disp('The solution is indeterminate');
0136 end
0137
0138 if unique
0139
0140 eu(2)=1;
0141 else
0142 fprintf(1,'Indeterminacy. %d loose endog errors.\n',k-r);
0143 disp('Generalized eigenvalues')
0144 disp(gev);
0145 disp(sort(abs(gev(:,2) ./ gev(:,1))));
0146
0147 end
0148
0149
0150 if isempty(Vhat1)
0151 unique=1;
0152 else
0153 unique=norm(Vhat1-V1*V1'*Vhat1)<realsmall*n;
0154 end
0155
0156 if unique
0157
0158
0159
0160 eu(2)=1;
0161 else
0162 fprintf(1,'Indeterminacy. %d loose endog errors.\n',size(Vhat1,2)-size(V1,2));
0163
0164
0165
0166
0167
0168
0169
0170 end
0171
0172
0173 PHI = (U1*(D11\V1')*Vhat1*Dhat11*Uhat1')';
0174
0175 tmat = [eye(n-nunstab) -PHI];
0176
0177 G0 = [ tmat * A ;
0178 zeros(m,n-m) eye(m)];
0179
0180 G1 = [ tmat * B ;
0181 zeros(m,n)];
0182
0183 G2 = [zeros(n-m,m);-eye(m)];
0184
0185 G0I=inv(G0);
0186
0187 H = Z*G0I;
0188
0189
0190
0191 if ~unique
0192 p = 1;
0193 m_1 = 0;
0194 m_2 = 1;
0195
0196 if m_1
0197 K1 = Q1*PPI-PHI*Q2*PPI;
0198 K2 = -V1*inv(D11)*U1'*Q2*PSI;
0199 AA = K1*K2;
0200 BB = K1*V2;
0201 M1 = inv(BB'*BB)*BB'*(-AA);
0202 else
0203 M1 = zeros(k-r,l);
0204 end
0205 if m_2
0206 M2 = ones(k-r,p);
0207 else
0208 M2 = zeros(k-r,p);
0209 end
0210 end
0211
0212
0213 S0 = real(H*[Q1 - PHI * Q2 ; inv(A22-B22)*Q2 ]*C);
0214 S1 = real(H*G1*Z');
0215 Sy = real(H*G2);
0216 Sf = real(B22\A22);
0217 Sz = real(B22\(Q2*PSI));
0218
0219 if unique
0220 S2 = H * [ Q1 - PHI * Q2 ; zeros(m,n) ] * PSI;
0221 S2 = real(S2);
0222
0223 S3 = zeros(0);
0224 else
0225 S2 = H*[Q1-PHI*Q2 ; zeros(m,n)]*PSI + H*[(Q1*PPI-PHI*Q2*PPI)*(-V1*inv(D11)*U1'*Q2*PSI+V2*M1); zeros(m,l)];
0226 S2 = real(S2);
0227
0228 S3 = H * [ (Q1*PPI - PHI*Q2*PPI)*V2*M2 ; zeros(m,p)];
0229 S3 = real(S3);
0230 end
0231
0232
0233
0234