From c51066d76d570708826c328a4b26af26222bca75 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Mon, 2 May 2011 11:12:39 +0200 Subject: [PATCH] Complete re-writing to make the function modular (the input being either J or H or Hess etc.). --- matlab/identification_checks.m | 196 +++++++++++---------------------- 1 file changed, 64 insertions(+), 132 deletions(-) diff --git a/matlab/identification_checks.m b/matlab/identification_checks.m index a3541207b..dd4778e8b 100644 --- a/matlab/identification_checks.m +++ b/matlab/identification_checks.m @@ -1,35 +1,22 @@ -function [McoH, McoJ, McoGP, PcoH, PcoJ, PcoGP, condH, condJ, condGP, eH, eJ, eGP, ind01, ind02, indnoH, indnoJ, ixnoH, ixnoJ] = identification_checks(H, JJ, gp) -% function [McoH, McoJ, McoGP, PcoH, PcoJ, PcoGP, condH, condJ, condGP, eH, -% eJ, eGP, ind01, ind02, indnoH, indnoJ, ixnoH, ixnoJ] = identification_checks(H, JJ, gp) +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) % checks for identification % % INPUTS -% o H [matrix] [(entries in st.sp. model solutio) x nparams] -% derivatives of model solution w.r.t. parameters and shocks -% o JJ [matrix] [moments x nparams] -% derivatives of moments w.r.t. parameters and shocks -% o gp [matrix] [jacobian_entries x nparams] -% derivatives of jacobian (i.e. LRE model) w.r.t. parameters and shocks +% 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 McoH [array] multicollinearity coefficients in the model solution -% o McoJ [array] multicollinearity coefficients in the moments -% o McoGP [array] multicollinearity coefficients in the LRE model -% o PcoH [matrix] pairwise correlations in the model solution -% o PcoJ [matrix] pairwise correlations in the moments -% o PcoGP [matrix] pairwise correlations in the LRE model -% o condH condition number of H -% o condJ condition number of JJ -% o condGP condition number of gp -% o eH eigevectors of H -% o eJ eigevectors of JJ -% o eGP eigevectors of gp -% o ind01 [array] binary indicator for non-zero columns of H -% o ind02 [array] binary indicator for non-zero columns of JJ -% o indnoH [matrix] index of non-identified params in H -% o indnoJ [matrix] index of non-identified params in JJ -% o ixnoH number of rows in ind01 -% o ixnoJ number of rows in ind02 +% 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 @@ -54,145 +41,90 @@ function [McoH, McoJ, McoGP, PcoH, PcoJ, PcoGP, condH, condJ, condGP, eH, eJ, eG % My suggestion is to have the following steps for identification check in % dynare: -% 1. check rank of H, JJ, gp at theta -npar = size(H,2); -npar0 = size(gp,2); % shocks do not enter jacobian -indnoH = zeros(1,npar); +% 1. check rank of JJ at theta +npar = size(JJ,2); indnoJ = zeros(1,npar); -indnoLRE = zeros(1,npar0); -% H matrix -ind1 = find(vnorm(H)>=eps); % take non-zero columns -H1 = H(:,ind1); -[eu,e2,e1] = svd( H1, 0 ); -eH = zeros(npar,npar); -% eH(ind1,:) = e1; -eH(ind1,length(find(vnorm(H)=eps); % take non-zero columns -JJ1 = JJ(:,ind2); +ind1 = find(vnorm(JJ)>=eps); % take non-zero columns +JJ1 = JJ(:,ind1); [eu,ee2,ee1] = svd( JJ1, 0 ); -eJ = zeros(npar,npar); -eJ(ind2,length(find(vnorm(JJ)=eps); % take non-zero columns -gp1 = gp(:,ind3); -covgp = gp1'*gp1; -sdgp = sqrt(diag(covgp)); -sdgp = sdgp*sdgp'; -[eu,ex2,ex1] = svd(gp1, 0 ); -eGP = zeros(npar0,npar0); -eGP(ind3,length(find(vnorm(gp) 1.e-3 )'; - 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(' ') +else + deltaJ = sqrt(diag(JJ)); + tildaJ = JJ./((deltaJ)*(deltaJ')); + McoJ(:,1)=(1-1./diag(inv(tildaJ))); + rhoM=sqrt(1-McoJ); +% PcoJ=inv(tildaJ); + PcoJ=inv(JJ); + sd=sqrt(diag(PcoJ)); + PcoJ = PcoJ./((sd)*(sd')); end -ixnoH=ixno; -ixno = 0; + +ixnoJ = 0; if rankJ 1.e-3)'; + ixnoJ = ixnoJ + 1; + indnoJ(ixnoJ,ind1) = (abs(ee1(:,ee0(j))) > 1.e-3)'; 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 % be checked using -PcoH = NaN(npar,npar); +jweak=zeros(1,npar); +jweak_pair=zeros(npar,npar); + +if hess_flag==0, PcoJ = NaN(npar,npar); -PcoGP = NaN(npar0,npar0); -for ii = 1:size(H1,2); - PcoH(ind1(ii),ind1(ii)) = 1; - for jj = ii+1:size(H1,2); - PcoH(ind1(ii),ind1(jj)) = [cosn([H1(:,ii),H1(:,jj)])]; - PcoH(ind1(jj),ind1(ii)) = PcoH(ind1(ii),ind1(jj)); - end -end for ii = 1:size(JJ1,2); - PcoJ(ind2(ii),ind2(ii)) = 1; + PcoJ(ind1(ii),ind1(ii)) = 1; for jj = ii+1:size(JJ1,2); - PcoJ(ind2(ii),ind2(jj)) = [cosn([JJ1(:,ii),JJ1(:,jj)])]; - PcoJ(ind2(jj),ind2(ii)) = PcoJ(ind2(ii),ind2(jj)); + PcoJ(ind1(ii),ind1(jj)) = cosn([JJ1(:,ii),JJ1(:,jj)]); + PcoJ(ind1(jj),ind1(ii)) = PcoJ(ind1(ii),ind1(jj)); end end -for ii = 1:size(gp1,2); - PcoGP(ind3(ii),ind3(ii)) = 1; - for jj = ii+1:size(gp1,2); - PcoGP(ind3(ii),ind3(jj)) = [cosn([gp1(:,ii),gp1(:,jj)])]; - PcoGP(ind3(jj),ind3(ii)) = PcoGP(ind3(ii),ind3(jj)); +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 end end +end - - - +jweak_pair=dyn_vech(jweak_pair)';