identification_checks.m 4.26 KB
Newer Older
1
2
function [condJ, ind0, indnoJ, ixnoJ, McoJ, PcoJ, jweak, jweak_pair] = identification_checks(JJ, hess_flag)
% function [condJ, ind0, indnoJ, ixnoJ, McoJ, PcoJ, jweak, jweak_pair] = identification_checks(JJ, hess_flag)
3
4
5
% checks for identification
%
% INPUTS
6
7
8
9
%    o JJ               [matrix] [output x nparams] IF hess_flag==0
%                                 derivatives of output w.r.t. parameters and shocks
%    o JJ               [matrix] [nparams x nparams] IF hess_flag==1
%                                 information matrix
10
11
%    
% OUTPUTS
12
13
14
15
16
17
18
19
%    o cond             condition number of JJ
%    o ind0             [array] binary indicator for non-zero columns of H
%    o indnoJ           [matrix] index of non-identified params 
%    o ixnoJ            number of rows in indnoJ
%    o Mco              [array] multicollinearity coefficients
%    o Pco              [matrix] pairwise correlations 
%    o jweak            [binary array] gives 1 if the  parameter has Mco=1(with tolerance 1.e-10)
%    o jweak_pair       [binary matrix] gives 1 if a couple parameters has Pco=1(with tolerance 1.e-10)
20
21
22
%    
% SPECIAL REQUIREMENTS
%    None
23

24
% Copyright (C) 2008-2011 Dynare Team
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
%
% 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/>.

% My suggestion is to have the following steps for identification check in
% dynare:

44
45
% 1. check rank of JJ at theta
npar = size(JJ,2);
46
indnoJ = zeros(1,npar);
47

48
49
50
51
52
if size(JJ,1)>1,
    ind1 = find(vnorm(JJ)>=eps); % take non-zero columns
else
    ind1 = find(abs(JJ)>=eps); % take non-zero columns
end    
53
JJ1 = JJ(:,ind1);
54
[eu,ee2,ee1] = svd( JJ1, 0 );
55
condJ= cond(JJ1);
56
rankJ = rank(JJ./norm(JJ),1.e-10);
57
rankJJ = rankJ;
58
59
60
% if hess_flag==0,
%     rankJJ = rank(JJ'*JJ);
% end   
61
62
63
64
65
66
67
68
69

ind0 = zeros(1,npar);
ind0(ind1) = 1;

if hess_flag==0,
    % find near linear dependence problems:
    McoJ = NaN(npar,1);
    for ii = 1:size(JJ1,2);
        McoJ(ind1(ii),:) = cosn([JJ1(:,ii),JJ1(:,find([1:1:size(JJ1,2)]~=ii))]);
70
    end
71
else
72
73
74
    deltaJ = sqrt(diag(JJ(ind1,ind1)));
    tildaJ = JJ(ind1,ind1)./((deltaJ)*(deltaJ'));
    McoJ(ind1,1)=(1-1./diag(inv(tildaJ)));
75
76
    rhoM=sqrt(1-McoJ);
%     PcoJ=inv(tildaJ);
77
78
    PcoJ=NaN(npar,npar);
    PcoJ(ind1,ind1)=inv(JJ(ind1,ind1));
79
    sd=sqrt(diag(PcoJ));
80
    PcoJ = abs(PcoJ./((sd)*(sd')));
81
82
end

83
84

ixnoJ = 0;
85
if rankJ<npar || rankJJ<npar || min(1-McoJ)<1.e-10
86
87
88
    %         - find out which parameters are involved
    %   disp('Some parameters are NOT identified by the moments included in J')
    %   disp(' ')
89
    if length(ind1)<npar,
90
        % parameters with zero column in JJ
91
92
        ixnoJ = ixnoJ + 1;
        indnoJ(ixnoJ,:) = (~ismember([1:npar],ind1));
93
    end
94
    ee0 = [rankJJ+1:length(ind1)];
95
    for j=1:length(ee0),
96
        % linearely dependent parameters in JJ
97
98
        ixnoJ = ixnoJ + 1;
        indnoJ(ixnoJ,ind1) = (abs(ee1(:,ee0(j))) > 1.e-3)';
99
100
    end
else  %rank(J)==length(theta) =>
101
      %         disp('All parameters are identified at theta by the moments included in J')
102
103
104
105
106
107
end

% here there is no exact linear dependence, but there are several
%     near-dependencies, mostly due to strong pairwise colliniearities, which can
%     be checked using

108
109
110
111
jweak=zeros(1,npar);
jweak_pair=zeros(npar,npar);

if hess_flag==0,
112
113
114
PcoJ = NaN(npar,npar);

for ii = 1:size(JJ1,2);
115
    PcoJ(ind1(ii),ind1(ii)) = 1;
116
    for jj = ii+1:size(JJ1,2);
117
118
        PcoJ(ind1(ii),ind1(jj)) = cosn([JJ1(:,ii),JJ1(:,jj)]);
        PcoJ(ind1(jj),ind1(ii)) = PcoJ(ind1(ii),ind1(jj));
119
120
121
    end
end

122
123
124
125
126
127
128
129
for j=1:npar,
    if McoJ(j)>(1-1.e-10), 
        jweak(j)=1;
        [ipair, jpair] = find(PcoJ(j,j+1:end)>(1-1.e-10));
        for jx=1:length(jpair),
            jweak_pair(j, jpair(jx)+j)=1;
            jweak_pair(jpair(jx)+j, j)=1;
        end
130
    end
131
end
132
end
133

134
jweak_pair=dyn_vech(jweak_pair)';
135