improved identification checks, using svd in place of eig.
parent
70400e5e84
commit
1fc3d503ee
|
@ -1,4 +1,4 @@
|
|||
function [McoH, McoJ, McoGP, PcoH, PcoJ, PcoGP, condH, condJ, condGP, eH, eJ, eGP, ind01, ind02, indnoH, indnoJ] = identification_checks(H,JJ, gp, bayestopt_)
|
||||
function [McoH, McoJ, McoGP, PcoH, PcoJ, PcoGP, condH, condJ, condGP, eH, eJ, eGP, ind01, ind02, indnoH, indnoJ, ixnoH, ixnoJ] = identification_checks(H,JJ, gp, bayestopt_)
|
||||
|
||||
% Copyright (C) 2008 Dynare Team
|
||||
%
|
||||
|
@ -23,36 +23,42 @@ function [McoH, McoJ, McoGP, PcoH, PcoJ, PcoGP, condH, condJ, condGP, eH, eJ, eG
|
|||
% 1. check rank of H at theta
|
||||
npar = size(H,2);
|
||||
npar0 = size(gp,2);
|
||||
indnoH = {};
|
||||
indnoJ = {};
|
||||
indnoLRE = {};
|
||||
ind1 = find(vnorm(H)~=0);
|
||||
indnoH = zeros(1,npar);
|
||||
indnoJ = zeros(1,npar);
|
||||
indnoLRE = zeros(1,npar0);
|
||||
ind1 = find(vnorm(H)>=eps);
|
||||
H1 = H(:,ind1);
|
||||
covH = H1'*H1;
|
||||
sdH = sqrt(diag(covH));
|
||||
sdH = sdH*sdH';
|
||||
[e1,e2] = eig( (H1'*H1)./sdH );
|
||||
% [e1,e2] = eig( (H1'*H1)./sdH );
|
||||
[eu,e2,e1] = svd( H1, 0 );
|
||||
eH = zeros(npar,npar);
|
||||
% eH(ind1,:) = e1;
|
||||
eH(ind1,length(find(vnorm(H)==0))+1:end) = e1;
|
||||
eH(find(vnorm(H)==0),1:length(find(vnorm(H)==0)))=eye(length(find(vnorm(H)==0)));
|
||||
condH = cond(H1);
|
||||
% condH = cond(H1'*H1);
|
||||
condHH = cond(H1'*H1);
|
||||
rankH = rank(H);
|
||||
rankHH = rank(H'*H);
|
||||
|
||||
ind2 = find(vnorm(JJ)~=0);
|
||||
ind2 = find(vnorm(JJ)>=eps);
|
||||
JJ1 = JJ(:,ind2);
|
||||
covJJ = JJ1'*JJ1;
|
||||
sdJJ = sqrt(diag(covJJ));
|
||||
sdJJ = sdJJ*sdJJ';
|
||||
[ee1,ee2] = eig( (JJ1'*JJ1)./sdJJ );
|
||||
% [ee1,ee2] = eig( (JJ1'*JJ1)./sdJJ );
|
||||
[eu,ee2,ee1] = svd( JJ1, 0 );
|
||||
% eJ = NaN(npar,length(ind2));
|
||||
eJ = zeros(npar,npar);
|
||||
eJ(ind2,length(find(vnorm(JJ)==0))+1:end) = ee1;
|
||||
eJ(find(vnorm(JJ)==0),1:length(find(vnorm(JJ)==0)))=eye(length(find(vnorm(JJ)==0)));
|
||||
% condJ = cond(JJ1'*JJ1);
|
||||
condJJ = cond(JJ1'*JJ1);
|
||||
condJ = cond(JJ1);
|
||||
rankJJ = rank(JJ'*JJ);
|
||||
rankJ = rank(JJ);
|
||||
|
||||
ind3 = find(vnorm(gp)~=0);
|
||||
ind3 = find(vnorm(gp)>=eps);
|
||||
gp1 = gp(:,ind3);
|
||||
covgp = gp1'*gp1;
|
||||
sdgp = sqrt(diag(covgp));
|
||||
|
@ -65,54 +71,11 @@ eGP(find(vnorm(gp)==0),1:length(find(vnorm(gp)==0)))=eye(length(find(vnorm(gp)==
|
|||
% condJ = cond(JJ1'*JJ1);
|
||||
condGP = cond(gp1);
|
||||
|
||||
if rank(H)<npar
|
||||
ixno = 0;
|
||||
% - find out which parameters are involved,
|
||||
% using something like the vnorm and the eigenvalue decomposition of H;
|
||||
% disp('Some parameters are NOT identified in the model: H rank deficient')
|
||||
% disp(' ')
|
||||
if length(ind1)<npar,
|
||||
ixno = ixno + 1;
|
||||
indnoH(ixno) = {find(~ismember([1:npar],ind1))};
|
||||
% disp('Not identified params')
|
||||
% disp(bayestopt_.name(indnoH{1}))
|
||||
% disp(' ')
|
||||
end
|
||||
e0 = find(abs(diag(e2))<eps);
|
||||
for j=1:length(e0),
|
||||
ixno = ixno + 1;
|
||||
indnoH(ixno) = {ind1(find(e1(:,e0(j))))};
|
||||
% disp('Perfectly collinear parameters')
|
||||
% disp(bayestopt_.name(indnoH{ixno}))
|
||||
% disp(' ')
|
||||
end
|
||||
else % rank(H)==length(theta), go to 2
|
||||
% 2. check rank of J
|
||||
% disp('All parameters are identified at theta in the model (rank of H)')
|
||||
% disp(' ')
|
||||
end
|
||||
|
||||
if rank(JJ)<npar
|
||||
ixno = 0;
|
||||
% - find out which parameters are involved
|
||||
% disp('Some parameters are NOT identified by the moments included in J')
|
||||
% disp(' ')
|
||||
if length(ind2)<npar,
|
||||
ixno = ixno + 1;
|
||||
indnoJ(ixno) = {find(~ismember([1:npar],ind2))};
|
||||
end
|
||||
ee0 = find(abs(diag(ee2))<eps);
|
||||
for j=1:length(ee0),
|
||||
ixno = ixno + 1;
|
||||
indnoJ(ixno) = {ind2(find(ee1(:,ee0(j))))};
|
||||
% disp('Perfectly collinear parameters in moments J')
|
||||
% disp(bayestopt_.name(indnoJ{ixno}))
|
||||
% disp(' ')
|
||||
end
|
||||
else %rank(J)==length(theta) =>
|
||||
% disp('All parameters are identified at theta by the moments included in J')
|
||||
end
|
||||
|
||||
ind01 = zeros(npar,1);
|
||||
ind02 = zeros(npar,1);
|
||||
ind01(ind1) = 1;
|
||||
ind02(ind2) = 1;
|
||||
|
||||
% rank(H1)==size(H1,2)
|
||||
% rank(JJ1)==size(JJ1,2)
|
||||
|
@ -135,6 +98,64 @@ end
|
|||
% format long % some are nearly 1
|
||||
% McoJ
|
||||
|
||||
ixno = 0;
|
||||
if rankH<npar | rankHH<npar | min(1-McoH)<1.e-10
|
||||
% - find out which parameters are involved,
|
||||
% using something like the vnorm and the eigenvalue decomposition of H;
|
||||
% disp('Some parameters are NOT identified in the model: H rank deficient')
|
||||
% disp(' ')
|
||||
if length(ind1)<npar,
|
||||
ixno = ixno + 1;
|
||||
% indnoH(ixno) = {find(~ismember([1:npar],ind1))};
|
||||
indnoH(ixno,:) = (~ismember([1:npar],ind1));
|
||||
% disp('Not identified params')
|
||||
% disp(bayestopt_.name(indnoH{1}))
|
||||
% disp(' ')
|
||||
end
|
||||
e0 = [rankHH+1:length(ind1)];
|
||||
for j=1:length(e0),
|
||||
ixno = ixno + 1;
|
||||
% indnoH(ixno) = {ind1(find(abs(e1(:,e0(j)))) > 1.e-6 )};
|
||||
indnoH(ixno,:) = (abs(e1(:,e0(j))) > 1.e-6 )';
|
||||
% disp('Perfectly collinear parameters')
|
||||
% disp(bayestopt_.name(indnoH{ixno}))
|
||||
% disp(' ')
|
||||
% ind01(indnoH{ixno})=0;
|
||||
end
|
||||
else % rank(H)==length(theta), go to 2
|
||||
% 2. check rank of J
|
||||
% disp('All parameters are identified at theta in the model (rank of H)')
|
||||
% disp(' ')
|
||||
end
|
||||
ixnoH=ixno;
|
||||
|
||||
ixno = 0;
|
||||
if rankJ<npar | rankJJ<npar | min(1-McoJ)<1.e-10
|
||||
% - find out which parameters are involved
|
||||
% disp('Some parameters are NOT identified by the moments included in J')
|
||||
% disp(' ')
|
||||
if length(ind2)<npar,
|
||||
ixno = ixno + 1;
|
||||
% indnoJ(ixno) = {find(~ismember([1:npar],ind2))};
|
||||
indnoJ(ixno,:) = (~ismember([1:npar],ind2));
|
||||
end
|
||||
ee0 = [rankJJ+1:length(ind2)];
|
||||
if isempty(ee0),
|
||||
cccc=0';
|
||||
end
|
||||
for j=1:length(ee0),
|
||||
ixno = ixno + 1;
|
||||
% indnoJ(ixno) = {ind2( find(abs(ee1(:,ee0(j))) > 1.e-6) )};
|
||||
indnoJ(ixno,:) = (abs(ee1(:,ee0(j))) > 1.e-6)';
|
||||
% disp('Perfectly collinear parameters in moments J')
|
||||
% disp(bayestopt_.name(indnoJ{ixno}))
|
||||
% disp(' ')
|
||||
% ind02(indnoJ{ixno})=0;
|
||||
end
|
||||
else %rank(J)==length(theta) =>
|
||||
% disp('All parameters are identified at theta by the moments included in J')
|
||||
end
|
||||
ixnoJ=ixno;
|
||||
|
||||
% here there is no exact linear dependence, but there are several
|
||||
% near-dependencies, mostly due to strong pairwise colliniearities, which can
|
||||
|
@ -168,10 +189,10 @@ for ii = 1:size(gp1,2);
|
|||
end
|
||||
|
||||
|
||||
ind01 = zeros(npar,1);
|
||||
ind02 = zeros(npar,1);
|
||||
ind01(ind1) = 1;
|
||||
ind02(ind2) = 1;
|
||||
% ind01 = zeros(npar,1);
|
||||
% ind02 = zeros(npar,1);
|
||||
% ind01(ind1) = 1;
|
||||
% ind02(ind2) = 1;
|
||||
|
||||
|
||||
|
||||
|
|
Loading…
Reference in New Issue