Home > matlab > identification_checks.m

identification_checks

PURPOSE ^

function [condJ, ind0, indnoJ, ixnoJ, McoJ, PcoJ, jweak, jweak_pair] = identification_checks(JJ, hess_flag)

SYNOPSIS ^

function [condJ, ind0, indnoJ, ixnoJ, McoJ, PcoJ, jweak, jweak_pair] = identification_checks(JJ, hess_flag)

DESCRIPTION ^

 function [condJ, ind0, indnoJ, ixnoJ, McoJ, PcoJ, jweak, jweak_pair] = identification_checks(JJ, hess_flag)
 checks for identification

 INPUTS
    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
    
 OUTPUTS
    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)
    
 SPECIAL REQUIREMENTS
    None

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [condJ, ind0, indnoJ, ixnoJ, McoJ, PcoJ, jweak, jweak_pair] = identification_checks(JJ, hess_flag)
0002 % function [condJ, ind0, indnoJ, ixnoJ, McoJ, PcoJ, jweak, jweak_pair] = identification_checks(JJ, hess_flag)
0003 % checks for identification
0004 %
0005 % INPUTS
0006 %    o JJ               [matrix] [output x nparams] IF hess_flag==0
0007 %                                 derivatives of output w.r.t. parameters and shocks
0008 %    o JJ               [matrix] [nparams x nparams] IF hess_flag==1
0009 %                                 information matrix
0010 %
0011 % OUTPUTS
0012 %    o cond             condition number of JJ
0013 %    o ind0             [array] binary indicator for non-zero columns of H
0014 %    o indnoJ           [matrix] index of non-identified params
0015 %    o ixnoJ            number of rows in indnoJ
0016 %    o Mco              [array] multicollinearity coefficients
0017 %    o Pco              [matrix] pairwise correlations
0018 %    o jweak            [binary array] gives 1 if the  parameter has Mco=1(with tolerance 1.e-10)
0019 %    o jweak_pair       [binary matrix] gives 1 if a couple parameters has Pco=1(with tolerance 1.e-10)
0020 %
0021 % SPECIAL REQUIREMENTS
0022 %    None
0023 
0024 % Copyright (C) 2008-2011 Dynare Team
0025 %
0026 % This file is part of Dynare.
0027 %
0028 % Dynare is free software: you can redistribute it and/or modify
0029 % it under the terms of the GNU General Public License as published by
0030 % the Free Software Foundation, either version 3 of the License, or
0031 % (at your option) any later version.
0032 %
0033 % Dynare is distributed in the hope that it will be useful,
0034 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0035 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0036 % GNU General Public License for more details.
0037 %
0038 % You should have received a copy of the GNU General Public License
0039 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0040 
0041 % My suggestion is to have the following steps for identification check in
0042 % dynare:
0043 
0044 % 1. check rank of JJ at theta
0045 npar = size(JJ,2);
0046 indnoJ = zeros(1,npar);
0047 
0048 ind1 = find(vnorm(JJ)>=eps); % take non-zero columns
0049 JJ1 = JJ(:,ind1);
0050 [eu,ee2,ee1] = svd( JJ1, 0 );
0051 condJ= cond(JJ1);
0052 rankJ = rank(JJ./norm(JJ),1.e-10);
0053 rankJJ = rankJ;
0054 % if hess_flag==0,
0055 %     rankJJ = rank(JJ'*JJ);
0056 % end
0057 
0058 ind0 = zeros(1,npar);
0059 ind0(ind1) = 1;
0060 
0061 if hess_flag==0,
0062     % find near linear dependence problems:
0063     McoJ = NaN(npar,1);
0064     for ii = 1:size(JJ1,2);
0065         McoJ(ind1(ii),:) = cosn([JJ1(:,ii),JJ1(:,find([1:1:size(JJ1,2)]~=ii))]);
0066     end
0067 else
0068     deltaJ = sqrt(diag(JJ(ind1,ind1)));
0069     tildaJ = JJ(ind1,ind1)./((deltaJ)*(deltaJ'));
0070     McoJ(ind1,1)=(1-1./diag(inv(tildaJ)));
0071     rhoM=sqrt(1-McoJ);
0072 %     PcoJ=inv(tildaJ);
0073     PcoJ=NaN(npar,npar);
0074     PcoJ(ind1,ind1)=inv(JJ(ind1,ind1));
0075     sd=sqrt(diag(PcoJ));
0076     PcoJ = abs(PcoJ./((sd)*(sd')));
0077 end
0078 
0079 
0080 ixnoJ = 0;
0081 if rankJ<npar || rankJJ<npar || min(1-McoJ)<1.e-10
0082     %         - find out which parameters are involved
0083     %   disp('Some parameters are NOT identified by the moments included in J')
0084     %   disp(' ')
0085     if length(ind1)<npar,
0086         % parameters with zero column in JJ
0087         ixnoJ = ixnoJ + 1;
0088         indnoJ(ixnoJ,:) = (~ismember([1:npar],ind1));
0089     end
0090     ee0 = [rankJJ+1:length(ind1)];
0091     for j=1:length(ee0),
0092         % linearely dependent parameters in JJ
0093         ixnoJ = ixnoJ + 1;
0094         indnoJ(ixnoJ,ind1) = (abs(ee1(:,ee0(j))) > 1.e-3)';
0095     end
0096 else  %rank(J)==length(theta) =>
0097       %         disp('All parameters are identified at theta by the moments included in J')
0098 end
0099 
0100 % here there is no exact linear dependence, but there are several
0101 %     near-dependencies, mostly due to strong pairwise colliniearities, which can
0102 %     be checked using
0103 
0104 jweak=zeros(1,npar);
0105 jweak_pair=zeros(npar,npar);
0106 
0107 if hess_flag==0,
0108 PcoJ = NaN(npar,npar);
0109 
0110 for ii = 1:size(JJ1,2);
0111     PcoJ(ind1(ii),ind1(ii)) = 1;
0112     for jj = ii+1:size(JJ1,2);
0113         PcoJ(ind1(ii),ind1(jj)) = cosn([JJ1(:,ii),JJ1(:,jj)]);
0114         PcoJ(ind1(jj),ind1(ii)) = PcoJ(ind1(ii),ind1(jj));
0115     end
0116 end
0117 
0118 for j=1:npar,
0119     if McoJ(j)>(1-1.e-10), 
0120         jweak(j)=1;
0121         [ipair, jpair] = find(PcoJ(j,j+1:end)>(1-1.e-10));
0122         for jx=1:length(jpair),
0123             jweak_pair(j, jpair(jx)+j)=1;
0124             jweak_pair(jpair(jx)+j, j)=1;
0125         end
0126     end
0127 end
0128 end
0129 
0130 jweak_pair=dyn_vech(jweak_pair)';
0131

Generated on Mon 21-May-2012 02:42:43 by m2html © 2005