stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge remote-tracking branch 'remotes/ratto/master'

time-shift
Sébastien Villemot 2011-05-02 11:35:17 +02:00
parent 19beebb28a cbfc017356
commit a0062605be
7 changed files with 901 additions and 631 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
matlab/cosn.m
View File

@ -27,14 +27,7 @@ function [co, b, yhat] = cosn(H);
y = H(:,1);
X = H(:,2:end);
% y = H(:,1);
% X = H(:,2:end);
if size(X,2)==1;
b=1;
else
b=(X\y);
end
b=(X\y);
yhat = X*b;
if rank(yhat),
co = abs(y'*yhat/sqrt((y'*y)*(yhat'*yhat)));

290
matlab/disp_identification.m
View File

@ -1,4 +1,4 @@
function disp_identification(pdraws, idemodel, idemoments, name, advanced)
function disp_identification(pdraws, idemodel, idemoments, name)
% Copyright (C) 2008-2011 Dynare Team
%
@ -19,69 +19,66 @@ function disp_identification(pdraws, idemodel, idemoments, name, advanced)
global options_
if nargin<5 || isempty(advanced),
advanced=0;
end
[SampleSize, npar] = size(pdraws);
jok = 0;
jokP = 0;
jokJ = 0;
jokPJ = 0;
% jok = 0;
% jokP = 0;
% jokJ = 0;
% jokPJ = 0;
for j=1:npar,
% if any(idemodel.ind(j,:)==0),
% pno = 100*length(find(idemodel.ind(j,:)==0))/SampleSize;
% disp(['Parameter ',name{j},' is not identified in the model for ',num2str(pno),'% of MC runs!' ])
% disp(' ')
% end
% if any(idemoments.ind(j,:)==0),
% pno = 100*length(find(idemoments.ind(j,:)==0))/SampleSize;
% disp(['Parameter ',name{j},' is not identified by J moments for ',num2str(pno),'% of MC runs!' ])
% disp(' ')
% end
if any(idemodel.ind(j,:)==1),
iok = find(idemodel.ind(j,:)==1);
jok = jok+1;
kok(jok) = j;
mmin(jok,1) = min(idemodel.Mco(j,iok));
mmean(jok,1) = mean(idemodel.Mco(j,iok));
mmax(jok,1) = max(idemodel.Mco(j,iok));
[ipmax, jpmax] = find(abs(squeeze(idemodel.Pco(j,[1:j-1,j+1:end],iok)))>0.95);
if ~isempty(ipmax)
jokP = jokP+1;
kokP(jokP) = j;
ipmax(find(ipmax>=j))=ipmax(find(ipmax>=j))+1;
[N,X]=hist(ipmax,[1:npar]);
jpM(jokP)={find(N)};
NPM(jokP)={N(find(N))./SampleSize.*100};
pmeanM(jokP)={mean(squeeze(idemodel.Pco(j,find(N),iok))')};
pminM(jokP)={min(squeeze(idemodel.Pco(j,find(N),iok))')};
pmaxM(jokP)={max(squeeze(idemodel.Pco(j,find(N),iok))')};
end
end
if any(idemoments.ind(j,:)==1),
iok = find(idemoments.ind(j,:)==1);
jokJ = jokJ+1;
kokJ(jokJ) = j;
mminJ(jokJ,1) = min(idemoments.Mco(j,iok));
mmeanJ(jokJ,1) = mean(idemoments.Mco(j,iok));
mmaxJ(jokJ,1) = max(idemoments.Mco(j,iok));
[ipmax, jpmax] = find(abs(squeeze(idemoments.Pco(j,[1:j-1,j+1:end],iok)))>0.95);
if ~isempty(ipmax)
jokPJ = jokPJ+1;
kokPJ(jokPJ) = j;
ipmax(find(ipmax>=j))=ipmax(find(ipmax>=j))+1;
[N,X]=hist(ipmax,[1:npar]);
jpJ(jokPJ)={find(N)};
NPJ(jokPJ)={N(find(N))./SampleSize.*100};
pmeanJ(jokPJ)={mean(squeeze(idemoments.Pco(j,find(N),iok))')};
pminJ(jokPJ)={min(squeeze(idemoments.Pco(j,find(N),iok))')};
pmaxJ(jokPJ)={max(squeeze(idemoments.Pco(j,find(N),iok))')};
end
end
end
% for j=1:npar,
% % if any(idemodel.ind(j,:)==0),
% % pno = 100*length(find(idemodel.ind(j,:)==0))/SampleSize;
% % disp(['Parameter ',name{j},' is not identified in the model for ',num2str(pno),'% of MC runs!' ])
% % disp(' ')
% % end
% % if any(idemoments.ind(j,:)==0),
% % pno = 100*length(find(idemoments.ind(j,:)==0))/SampleSize;
% % disp(['Parameter ',name{j},' is not identified by J moments for ',num2str(pno),'% of MC runs!' ])
% % disp(' ')
% % end
% if any(idemodel.ind(j,:)==1),
% iok = find(idemodel.ind(j,:)==1);
% jok = jok+1;
% kok(jok) = j;
% mmin(jok,1) = min(idemodel.Mco(j,iok));
% mmean(jok,1) = mean(idemodel.Mco(j,iok));
% mmax(jok,1) = max(idemodel.Mco(j,iok));
% [ipmax, jpmax] = find(abs(squeeze(idemodel.Pco(j,[1:j-1,j+1:end],iok)))>0.95);
% if ~isempty(ipmax)
% jokP = jokP+1;
% kokP(jokP) = j;
% ipmax(find(ipmax>=j))=ipmax(find(ipmax>=j))+1;
% [N,X]=hist(ipmax,[1:npar]);
% jpM(jokP)={find(N)};
% NPM(jokP)={N(find(N))./SampleSize.*100};
% pmeanM(jokP)={mean(squeeze(idemodel.Pco(j,find(N),iok))')};
% pminM(jokP)={min(squeeze(idemodel.Pco(j,find(N),iok))')};
% pmaxM(jokP)={max(squeeze(idemodel.Pco(j,find(N),iok))')};
% end
% end
% if any(idemoments.ind(j,:)==1),
% iok = find(idemoments.ind(j,:)==1);
% jokJ = jokJ+1;
% kokJ(jokJ) = j;
% mminJ(jokJ,1) = min(idemoments.Mco(j,iok));
% mmeanJ(jokJ,1) = mean(idemoments.Mco(j,iok));
% mmaxJ(jokJ,1) = max(idemoments.Mco(j,iok));
% [ipmax, jpmax] = find(abs(squeeze(idemoments.Pco(j,[1:j-1,j+1:end],iok)))>0.95);
% if ~isempty(ipmax)
% jokPJ = jokPJ+1;
% kokPJ(jokPJ) = j;
% ipmax(find(ipmax>=j))=ipmax(find(ipmax>=j))+1;
% [N,X]=hist(ipmax,[1:npar]);
% jpJ(jokPJ)={find(N)};
% NPJ(jokPJ)={N(find(N))./SampleSize.*100};
% pmeanJ(jokPJ)={mean(squeeze(idemoments.Pco(j,find(N),iok))')};
% pminJ(jokPJ)={min(squeeze(idemoments.Pco(j,find(N),iok))')};
% pmaxJ(jokPJ)={max(squeeze(idemoments.Pco(j,find(N),iok))')};
% end
% end
% end
disp([' ']),
if any(idemodel.ino),
disp('WARNING !!!')
@ -90,76 +87,120 @@ if any(idemodel.ino),
else
disp(['The rank of H (model) is deficient!' ]),
end
end
for j=1:npar,
if any(idemodel.ind(j,:)==0),
pno = 100*length(find(idemodel.ind(j,:)==0))/SampleSize;
if SampleSize>1
disp([name{j},' is not identified in the model for ',num2str(pno),'% of MC runs!' ])
else
disp([name{j},' is not identified in the model!' ])
disp(' ')
for j=1:npar,
if any(idemodel.ind0(:,j)==0),
pno = 100*length(find(idemodel.ind0(:,j)==0))/SampleSize;
if SampleSize>1
disp([' ',name{j},' is not identified in the model for ',num2str(pno),'% of MC runs!' ])
else
disp([' ',name{j},' is not identified in the model!' ])
end
disp([' [dJ/d(',name{j},')=0 for all tau elements in the model solution!]' ])
end
end
iweak = length(find(idemodel.Mco(j,:)'>(1-1.e-10)));
if iweak,
% disp('WARNING !!!')
% disp(['Model derivatives of parameter ',name{j},' are multi-collinear (with tol = 1.e-10) for ',num2str(iweak/SampleSize*100),'% of MC runs!' ])
if SampleSize>1
disp([name{j},' is collinear w.r.t. all other params ',num2str(iweak/SampleSize*100),'% of MC runs!' ])
else
disp([name{j},' is collinear w.r.t. all other params!' ])
npairs=size(idemodel.jweak_pair,2);
jmap_pair=dyn_unvech(1:npairs);
jstore=[];
disp(' ')
for j=1:npairs,
iweak = length(find(idemodel.jweak_pair(:,j)));
if iweak,
[jx,jy]=find(jmap_pair==j);
jstore=[jstore jx(1) jy(1)];
if SampleSize > 1
disp([' [',name{jx(1)},',',name{jy(1)},'] are PAIRWISE collinear (with tol = 1.e-10) for ',num2str(length(iweak)/SampleSize*100),'% of MC runs!' ])
else
disp([' [',name{jx(1)},',',name{jy(1)},'] are PAIRWISE collinear (with tol = 1.e-10) !' ])
end
end
if npar>(j+1),
[ipair, jpair] = find(squeeze(idemodel.Pco(j,j+1:end,:))'>(1-1.e-10));
else
[ipair, jpair] = find(squeeze(idemodel.Pco(j,j+1:end,:))>(1-1.e-10));
end
if ~isempty(jpair),
for jx=j+1:npar,
ixp = find(jx==(jpair+j));
if ~isempty(ixp)
if SampleSize > 1,
disp([' [',name{j},',',name{jx},'] are PAIRWISE collinear (with tol = 1.e-10) for ',num2str(length(ixp)/SampleSize*100),'% of MC runs!' ])
else
disp([' [',name{j},',',name{jx},'] are PAIRWISE collinear (with tol = 1.e-10)!' ])
end
end
end
disp(' ')
for j=1:npar,
iweak = length(find(idemodel.jweak(:,j)));
if iweak && ~ismember(j,jstore),
% disp('WARNING !!!')
% disp(['Model derivatives of parameter ',name{j},' are multi-collinear (with tol = 1.e-10) for ',num2str(iweak/SampleSize*100),'% of MC runs!' ])
if SampleSize>1
disp([name{j},' is collinear w.r.t. all other params ',num2str(iweak/SampleSize*100),'% of MC runs!' ])
else
disp([name{j},' is collinear w.r.t. all other params!' ])
end
end
end
% if npar>(j+1),
% [ipair, jpair] = find(squeeze(idemodel.Pco(j,j+1:end,:))'>(1-1.e-10));
% else
% [ipair, jpair] = find(squeeze(idemodel.Pco(j,j+1:end,:))>(1-1.e-10));
% end
% if ~isempty(jpair),
% for jx=j+1:npar,
% ixp = find(jx==(jpair+j));
% if ~isempty(ixp)
% if SampleSize > 1,
% disp([' [',name{j},',',name{jx},'] are PAIRWISE collinear (with tol = 1.e-10) for ',num2str(length(ixp)/SampleSize*100),'% of MC runs!' ])
% else
% disp([' [',name{j},',',name{jx},'] are PAIRWISE collinear (with tol = 1.e-10)!' ])
% end
% end
% end
% end
end
if ~any(idemodel.ino) && ~any(any(idemodel.ind==0))
if ~any(idemodel.ino) && ~any(any(idemodel.ind0==0))
disp(['All parameters are identified in the model (rank of H).' ]),
disp(' ')
end
if any(idemoments.ino),
disp(' ')
disp('WARNING !!!')
if SampleSize > 1,
disp(['The rank of J (moments) is deficient for ', num2str(length(find(idemoments.ino))/SampleSize*100),'% of MC runs!' ]),
else
disp(['The rank of J (moments) is deficient!' ]),
end
end
if any(idemoments.ino),
% disp('WARNING !!!')
% disp(['The rank of J (moments) is deficient for ', num2str(length(find(idemoments.ino))/SampleSize*100),'% of MC runs!' ]),
indno=[];
for j=1:SampleSize, indno=[indno;idemoments.indno{j}]; end
freqno = mean(indno)*100;
ifreq=find(freqno);
% indno=[];
% for j=1:SampleSize, indno=[indno;idemoments.indno{j}]; end
% freqno = mean(indno)*100;
% ifreq=find(freqno);
% disp('MOMENT RANK FAILURE DUE TO COLLINEARITY OF PARAMETERS:');
disp(' ')
for j=1:npar,
if any(idemoments.ind(j,:)==0),
pno = 100*length(find(idemoments.ind(j,:)==0))/SampleSize;
if any(idemoments.ind0(:,j)==0),
pno = 100*length(find(idemoments.ind0(:,j)==0))/SampleSize;
if SampleSize > 1
disp([name{j},' is not identified by J moments for ',num2str(pno),'% of MC runs!' ])
disp([' ',name{j},' is not identified by J moments for ',num2str(pno),'% of MC runs!' ])
else
disp([name{j},' is not identified by J moments!' ])
disp([' ',name{j},' is not identified by J moments!' ])
end
disp([' [dJ/d(',name{j},')=0 for all J moments!]' ])
end
end
disp(' ')
npairs=size(idemoments.jweak_pair,2);
jmap_pair=dyn_unvech(1:npairs);
jstore=[];
for j=1:npairs,
iweak = length(find(idemoments.jweak_pair(:,j)));
if iweak,
[jx,jy]=find(jmap_pair==j);
jstore=[jstore jx(1) jy(1)]';
if SampleSize > 1
disp([' [',name{jx(1)},',',name{jy(1)},'] are PAIRWISE collinear (with tol = 1.e-10) for ',num2str(length(iweak)/SampleSize*100),'% of MC runs!' ])
else
disp([' [',name{jx(1)},',',name{jy(1)},'] are PAIRWISE collinear (with tol = 1.e-10) !' ])
end
end
iweak = length(find(idemoments.Mco(j,:)'>(1-1.e-10)));
if iweak,
end
disp(' ')
for j=1:npar,
iweak = length(find(idemoments.jweak(:,j)));
if iweak && ~ismember(j,jstore),
% disp('WARNING !!!')
% disp(['Moment derivatives of parameter ',name{j},' are multi-collinear (with tol = 1.e-10) for ',num2str(iweak/SampleSize*100),'% of MC runs!' ])
if SampleSize > 1,
@ -167,27 +208,30 @@ if any(idemoments.ino),
else
disp([name{j},' is collinear w.r.t. all other params!' ])
end
if npar>(j+1),
[ipair, jpair] = find(squeeze(idemoments.Pco(j,j+1:end,:))'>(1-1.e-10));
else
[ipair, jpair] = find(squeeze(idemoments.Pco(j,j+1:end,:))>(1-1.e-10));
end
if ~isempty(jpair),
for jx=j+1:npar,
ixp = find(jx==(jpair+j));
if ~isempty(ixp)
if SampleSize > 1
disp([' [',name{j},',',name{jx},'] are PAIRWISE collinear (with tol = 1.e-10) for ',num2str(length(ixp)/SampleSize*100),'% of MC runs!' ])
else
disp([' [',name{j},',',name{jx},'] are PAIRWISE collinear (with tol = 1.e-10) !' ])
end
end
end
end
end
end
% if npar>(j+1),
% [ipair, jpair] = find(squeeze(idemoments.Pco(j,j+1:end,:))'>(1-1.e-10));
% else
% [ipair, jpair] = find(squeeze(idemoments.Pco(j,j+1:end,:))>(1-1.e-10));
% end
% if ~isempty(jpair),
% for jx=j+1:npar,
% ixp = find(jx==(jpair+j));
% if ~isempty(ixp)
% if SampleSize > 1
% disp([' [',name{j},',',name{jx},'] are PAIRWISE collinear (with tol = 1.e-10) for ',num2str(length(ixp)/SampleSize*100),'% of MC runs!' ])
% else
% disp([' [',name{j},',',name{jx},'] are PAIRWISE collinear (with tol = 1.e-10) !' ])
% end
% end
% end
% end
% end
% end
end
if ~any(idemoments.ino) && ~any(any(idemoments.ind==0))
if ~any(idemoments.ino) && ~any(any(idemoments.ind0==0))
disp(' ')
disp(['All parameters are identified by J moments (rank of J)' ]),
disp(' ')
end

540
matlab/dynare_identification.m
View File

@ -64,7 +64,9 @@ else
prior_exist=1;
end
% options_ident.load_ident_files=1;
iload = options_ident.load_ident_files;
%options_ident.advanced=1;
advanced = options_ident.advanced;
nlags = options_ident.ar;
periods = options_ident.periods;
@ -138,8 +140,8 @@ else
name_tex = [cellstr(M_.exo_names_tex); cellstr(M_.param_names_tex)];
end
options_ident = set_default_option(options_ident,'max_ord_bruteforce',min([2,nparam-1]));
max_ord_bruteforce = options_ident.max_ord_bruteforce;
options_ident = set_default_option(options_ident,'max_dim_cova_group',min([2,nparam-1]));
max_dim_cova_group = options_ident.max_dim_cova_group;
MaxNumberOfBytes=options_.MaxNumberOfBytes;
@ -150,192 +152,158 @@ disp(' ')
if iload <=0,
iteration = 0;
loop_indx = 0;
file_index = 0;
run_index = 0;
[I,J]=find(M_.lead_lag_incidence');
if prior_exist,
if exist([fname_,'_mean.mat'],'file'),
% disp('Testing posterior mean')
load([fname_,'_mean'],'xparam1')
pmean = xparam1';
clear xparam1
end
if exist([fname_,'_mode.mat'],'file'),
% disp('Testing posterior mode')
load([fname_,'_mode'],'xparam1')
pmode = xparam1';
clear xparam1
end
disp('Testing prior mean')
params = set_prior(estim_params_,M_,options_)';
else
params = [sqrt(diag(M_.Sigma_e))', M_.params'];
end
[idehess_prior, idemoments_prior, idemodel_prior, idelre_prior, derivatives_info_prior] = ...
identification_analysis(params,indx,indexo,options_ident,data_info, prior_exist, name_tex,1);
idehess_prior.params=params;
% siH = idemodel_prior.siH;
% siJ = idemoments_prior.siJ;
% siLRE = idelre_prior.siLRE;
% normH = max(abs(siH)')';
% normJ = max(abs(siJ)')';
% normLRE = max(abs(siLRE)')';
save([IdentifDirectoryName '/' M_.fname '_identif.mat'], 'idehess_prior', 'idemoments_prior','idemodel_prior', 'idelre_prior')
disp_identification(params, idemodel_prior, idemoments_prior, name);
plot_identification(params,idemoments_prior,idehess_prior,idemodel_prior,idelre_prior,advanced,'Prior mean - ',name);
if exist('pmode','var'),
disp('Testing posterior mode')
[idehess_mode, idemoments_mode, idemodel_mode, idelre_mode, derivatives_info_mode] = ...
identification_analysis(pmode,indx,indexo,options_ident,data_info, prior_exist, name_tex,0);
idehess_mode.params=pmode;
save([IdentifDirectoryName '/' M_.fname '_identif.mat'], 'idehess_mode', 'idemoments_mode','idemodel_mode', 'idelre_mode','-append')
disp_identification(params, idemodel_mode, idemoments_mode, name);
plot_identification(params,idemoments_mode,idehess_mode,idemodel_mode,idelre_mode,advanced,'Posterior mode - ',name);
end
if exist('pmean','var'),
disp('Testing posterior mean')
[idehess_mean, idemoments_mean, idemodel_mean, idelre_mean, derivatives_info_mean] = ...
identification_analysis(pmean,indx,indexo,options_ident,data_info, prior_exist, name_tex,0);
idehess_mean.params=pmean;
save([IdentifDirectoryName '/' M_.fname '_identif.mat'], 'idehess_mean', 'idemoments_mean','idemodel_mean', 'idelre_mean','-append')
disp_identification(params, idemodel_mean, idemoments_mean, name);
plot_identification(params,idemoments_mean,idehess_mean,idemodel_mean,idelre_mean,advanced,'Posterior mean - ',name);
end
if SampleSize > 1,
disp(' ')
disp('Monte Carlo Testing')
h = waitbar(0,'Monte Carlo identification checks ...');
iteration = 0;
loop_indx = 0;
file_index = 0;
run_index = 0;
options_MC=options_ident;
options_MC.advanced=0;
else
iteration = 1;
pdraws = [];
end
[I,J]=find(M_.lead_lag_incidence');
while iteration < SampleSize,
loop_indx = loop_indx+1;
if prior_exist,
if SampleSize==1,
if exist([fname_,'_mean.mat'],'file'),
disp('Testing posterior mean')
load([fname_,'_mean'],'xparam1')
params = xparam1';
clear xparam1
elseif exist([fname_,'_mode.mat'],'file'),
disp('Testing posterior mode')
load([fname_,'_mode'],'xparam1')
params = xparam1';
clear xparam1
else
disp('Testing prior mean')
params = set_prior(estim_params_,M_,options_)';
end
else
if nargin==2,
params = pdraws0(iteration+1,:);
else
params = prior_draw();
end
end
set_all_parameters(params);
if nargin==2,
params = pdraws0(iteration+1,:);
else
params = [sqrt(diag(M_.Sigma_e))', M_.params'];
params = prior_draw();
end
[~, ideJ, ideH, ideGP, ~ , info] = ...
identification_analysis(params,indx,indexo,options_MC,data_info, prior_exist, name_tex,0);
if iteration==0,
MAX_tau = min(SampleSize,ceil(MaxNumberOfBytes/(size(ideH.siH,1)*nparam)/8));
stoH = zeros([size(ideH.siH,1),nparam,MAX_tau]);
stoJJ = zeros([size(ideJ.siJ,1),nparam,MAX_tau]);
stoLRE = zeros([size(ideGP.siLRE,1),np,MAX_tau]);
TAU = zeros(size(ideH.siH,1),SampleSize);
GAM = zeros(size(ideJ.siJ,1),SampleSize);
LRE = zeros(size(ideGP.siLRE,1),SampleSize);
pdraws = zeros(SampleSize,nparam);
idemoments.indJJ = ideJ.indJJ;
idemodel.indH = ideH.indH;
idelre.indLRE = ideGP.indLRE;
idemoments.ind0 = zeros(SampleSize,nparam);
idemodel.ind0 = zeros(SampleSize,nparam);
idelre.ind0 = zeros(SampleSize,np);
idemoments.jweak = zeros(SampleSize,nparam);
idemodel.jweak = zeros(SampleSize,nparam);
idelre.jweak = zeros(SampleSize,np);
idemoments.jweak_pair = zeros(SampleSize,nparam*(nparam+1)/2);
idemodel.jweak_pair = zeros(SampleSize,nparam*(nparam+1)/2);
idelre.jweak_pair = zeros(SampleSize,np*(np+1)/2);
idemoments.cond = zeros(SampleSize,1);
idemodel.cond = zeros(SampleSize,1);
idelre.cond = zeros(SampleSize,1);
idemoments.Mco = zeros(SampleSize,nparam);
idemodel.Mco = zeros(SampleSize,nparam);
idelre.Mco = zeros(SampleSize,np);
delete([IdentifDirectoryName '/' M_.fname '_identif_*.mat'])
end
[A,B,ys,info]=dynare_resolve;
if info(1)==0,
oo0=oo_;
tau=[oo_.dr.ys(oo_.dr.order_var); vec(A); dyn_vech(B*M_.Sigma_e*B')];
yy0=oo_.dr.ys(I);
[residual, g1 ] = feval([M_.fname,'_dynamic'],yy0, ...
oo_.exo_steady_state', M_.params, ...
oo_.dr.ys, 1);
iteration = iteration + 1;
run_index = run_index + 1;
if iteration,
[JJ, H, gam, gp, dA, dOm, dYss] = getJJ(A, B, M_,oo0,options_,0,indx,indexo,bayestopt_.mf2,nlags,useautocorr);
derivatives_info.DT=dA;
derivatives_info.DOm=dOm;
derivatives_info.DYss=dYss;
if iteration==1,
indJJ = (find(max(abs(JJ'))>1.e-8));
indH = (find(max(abs(H'))>1.e-8));
indLRE = (find(max(abs(gp'))>1.e-8));
TAU = zeros(length(indH),SampleSize);
GAM = zeros(length(indJJ),SampleSize);
LRE = zeros(length(indLRE),SampleSize);
MAX_tau = min(SampleSize,ceil(MaxNumberOfBytes/(length(indH)*nparam)/8));
MAX_gam = min(SampleSize,ceil(MaxNumberOfBytes/(length(indJJ)*nparam)/8));
stoH = zeros([length(indH),nparam,MAX_tau]);
stoJJ = zeros([length(indJJ),nparam,MAX_tau]);
delete([IdentifDirectoryName '/' M_.fname '_identif_*.mat'])
end
TAU(:,iteration)=tau(indH);
vg1 = [oo_.dr.ys(oo_.dr.order_var); vec(g1)];
LRE(:,iteration)=vg1(indLRE);
GAM(:,iteration)=gam(indJJ);
stoLRE(:,:,run_index) = gp(indLRE,:);
stoH(:,:,run_index) = H(indH,:);
stoJJ(:,:,run_index) = JJ(indJJ,:);
% use prior uncertainty
siJ = (JJ(indJJ,:));
siH = (H(indH,:));
siLRE = (gp(indLRE,:));
pdraws(iteration,:) = params;
normH = max(abs(stoH(:,:,run_index))')';
normJ = max(abs(stoJJ(:,:,run_index))')';
normLRE = max(abs(stoLRE(:,:,run_index))')';
[idemodel.Mco(:,iteration), idemoments.Mco(:,iteration), idelre.Mco(:,iteration), ...
idemodel.Pco(:,:,iteration), idemoments.Pco(:,:,iteration), idelre.Pco(:,:,iteration), ...
idemodel.cond(iteration), idemoments.cond(iteration), idelre.cond(iteration), ...
idemodel.ee(:,:,iteration), idemoments.ee(:,:,iteration), idelre.ee(:,:,iteration), ...
idemodel.ind(:,iteration), idemoments.ind(:,iteration), ...
idemodel.indno{iteration}, idemoments.indno{iteration}, ...
idemodel.ino(iteration), idemoments.ino(iteration)] = ...
identification_checks(H(indH,:)./normH(:,ones(nparam,1)),JJ(indJJ,:)./normJ(:,ones(nparam,1)), gp(indLRE,:)./normLRE(:,ones(size(gp,2),1)));
% identification_checks(H(indH,:),JJ(indJJ,:), gp(indLRE,:), bayestopt_);
indok = find(max(idemoments.indno{iteration},[],1)==0);
if iteration ==1,
ide_strength_J=NaN(1,nparam);
ide_strength_J_prior=NaN(1,nparam);
end
if iteration ==1 && advanced,
[pars, cosnJ] = ident_bruteforce(JJ(indJJ,:)./normJ(:,ones(nparam,1)),max_ord_bruteforce,options_.TeX,name_tex);
end
if iteration ==1 && ~isempty(indok),
normaliz = abs(params);
if prior_exist,
if ~isempty(estim_params_.var_exo),
normaliz1 = estim_params_.var_exo(:,7); % normalize with prior standard deviation
else
normaliz1=[];
end
if ~isempty(estim_params_.param_vals),
normaliz1 = [normaliz1; estim_params_.param_vals(:,7)]'; % normalize with prior standard deviation
end
% normaliz = max([normaliz; normaliz1]);
normaliz1(isinf(normaliz1)) = 1;
else
normaliz1 = ones(1,nparam);
end
replic = max([replic, length(indJJ)*3]);
try,
options_.irf = 0;
options_.noprint = 1;
options_.order = 1;
options_.periods = data_info.gend+100;
info = stoch_simul(options_.varobs);
datax=oo_.endo_simul(options_.varobs_id,100+1:end);
% datax=data;
derivatives_info.no_DLIK=1;
[fval,cost_flag,ys,trend_coeff,info,DLIK,AHess] = DsgeLikelihood(params',data_info.gend,datax,data_info.data_index,data_info.number_of_observations,data_info.no_more_missing_observations,derivatives_info);
cparam = inv(-AHess);
normaliz(indok) = sqrt(diag(cparam))';
cmm = siJ*((-AHess)\siJ');
flag_score=1;
catch,
cmm = simulated_moment_uncertainty(indJJ, periods, replic);
% Jinv=(siJ(:,indok)'*siJ(:,indok))\siJ(:,indok)';
% MIM=inv(Jinv*cmm*Jinv');
MIM=siJ(:,indok)'*(cmm\siJ(:,indok));
deltaM = sqrt(diag(MIM));
tildaM = MIM./((deltaM)*(deltaM'));
rhoM=sqrt(1-1./diag(inv(tildaM)));
deltaM = deltaM.*normaliz(indok)';
normaliz(indok) = sqrt(diag(inv(MIM)))';
flag_score=0;
end
ide_strength_J(indok) = (1./(normaliz(indok)'./abs(params(indok)')));
ide_strength_J_prior(indok) = (1./(normaliz(indok)'./normaliz1(indok)'));
ide_strength_J(params==0)=ide_strength_J_prior(params==0);
quant = siJ./repmat(sqrt(diag(cmm)),1,nparam);
siJnorm(iteration,:) = vnorm(quant).*normaliz;
% siJnorm(iteration,:) = vnorm(siJ(inok,:)).*normaliz;
quant=[];
inok = find((abs(TAU(:,iteration))<1.e-8));
isok = find((abs(TAU(:,iteration))));
quant(isok,:) = siH(isok,:)./repmat(TAU(isok,iteration),1,nparam);
quant(inok,:) = siH(inok,:)./repmat(mean(abs(TAU(:,iteration))),length(inok),nparam);
siHnorm(iteration,:) = vnorm(quant).*normaliz;
% siHnorm(iteration,:) = vnorm(siH./repmat(TAU(:,iteration),1,nparam)).*normaliz;
quant=[];
inok = find((abs(LRE(:,iteration))<1.e-8));
isok = find((abs(LRE(:,iteration))));
quant(isok,:) = siLRE(isok,:)./repmat(LRE(isok,iteration),1,np);
quant(inok,:) = siLRE(inok,:)./repmat(mean(abs(LRE(:,iteration))),length(inok),np);
siLREnorm(iteration,:) = vnorm(quant).*normaliz(offset+1:end);
% siLREnorm(iteration,:) = vnorm(siLRE./repmat(LRE(:,iteration),1,nparam-offset)).*normaliz(offset+1:end);
end,
if run_index==MAX_tau || iteration==SampleSize,
file_index = file_index + 1;
if run_index<MAX_tau,
stoH = stoH(:,:,1:run_index);
stoJJ = stoJJ(:,:,1:run_index);
stoLRE = stoLRE(:,:,1:run_index);
end
save([IdentifDirectoryName '/' M_.fname '_identif_' int2str(file_index) '.mat'], 'stoH', 'stoJJ', 'stoLRE')
run_index = 0;
TAU(:,iteration)=ideH.TAU;
LRE(:,iteration)=ideGP.LRE;
GAM(:,iteration)=ideJ.GAM;
idemoments.cond(iteration,1)=ideJ.cond;
idemodel.cond(iteration,1)=ideH.cond;
idelre.cond(iteration,1)=ideGP.cond;
idemoments.ino(iteration,1)=ideJ.ino;
idemodel.ino(iteration,1)=ideH.ino;
idelre.ino(iteration,1)=ideGP.ino;
idemoments.ind0(iteration,:)=ideJ.ind0;
idemodel.ind0(iteration,:)=ideH.ind0;
idelre.ind0(iteration,:)=ideGP.ind0;
idemoments.jweak(iteration,:)=ideJ.jweak;
idemodel.jweak(iteration,:)=ideH.jweak;
idelre.jweak(iteration,:)=ideGP.jweak;
idemoments.jweak_pair(iteration,:)=ideJ.jweak_pair;
idemodel.jweak_pair(iteration,:)=ideH.jweak_pair;
idelre.jweak_pair(iteration,:)=ideGP.jweak_pair;
idemoments.Mco(iteration,:)=ideJ.Mco;
idemodel.Mco(iteration,:)=ideH.Mco;
idelre.Mco(iteration,:)=ideGP.Mco;
stoLRE(:,:,run_index) = ideGP.siLRE;
stoH(:,:,run_index) = ideH.siH;
stoJJ(:,:,run_index) = ideJ.siJ;
pdraws(iteration,:) = params;
if run_index==MAX_tau || iteration==SampleSize,
file_index = file_index + 1;
if run_index<MAX_tau,
stoH = stoH(:,:,1:run_index);
stoJJ = stoJJ(:,:,1:run_index);
stoLRE = stoLRE(:,:,1:run_index);
end
save([IdentifDirectoryName '/' M_.fname '_identif_' int2str(file_index) '.mat'], 'stoH', 'stoJJ', 'stoLRE')
run_index = 0;
stoH = zeros(size(stoH));
stoJJ = zeros(size(stoJJ));
stoLRE = zeros(size(stoLRE));
if SampleSize > 1,
waitbar(iteration/SampleSize,h,['MC Identification checks ',int2str(iteration),'/',int2str(SampleSize)])
end
end
if SampleSize > 1,
waitbar(iteration/SampleSize,h,['MC identification checks ',int2str(iteration),'/',int2str(SampleSize)])
end
end
end
@ -356,13 +324,19 @@ if iload <=0,
end
end
idemoments.siJnorm = siJnorm;
idemodel.siHnorm = siHnorm;
idelre.siLREnorm = siLREnorm;
save([IdentifDirectoryName '/' M_.fname '_identif.mat'], 'pdraws', 'idemodel', 'idemoments', 'idelre', ... %'indJJ', 'indH', 'indLRE', ...
'TAU', 'GAM', 'LRE','-append')
else
siJnorm = idemoments_prior.siJnorm;
siHnorm = idemodel_prior.siHnorm;
siLREnorm = idelre_prior.siLREnorm;
end
save([IdentifDirectoryName '/' M_.fname '_identif.mat'], 'pdraws', 'idemodel', 'idemoments', 'idelre', 'indJJ', 'indH', 'indLRE', ...
'siJnorm', 'siHnorm', 'siLREnorm', 'TAU', 'GAM', 'LRE')
else
load([IdentifDirectoryName '/' M_.fname '_identif'], 'pdraws', 'idemodel', 'idemoments', 'idelre', 'indJJ', 'indH', 'indLRE', ...
'siJnorm', 'siHnorm', 'siLREnorm', 'TAU', 'GAM', 'LRE')
load([IdentifDirectoryName '/' M_.fname '_identif'])
identFiles = dir([IdentifDirectoryName '/' M_.fname '_identif_*']);
options_ident.prior_mc=size(pdraws,1);
SampleSize = options_ident.prior_mc;
@ -383,202 +357,36 @@ if nargout>3 && iload,
end
end
disp_identification(pdraws, idemodel, idemoments, name, advanced)
figure('Name','Identification in the moments'),
subplot(211)
mmm = (ide_strength_J);
[ss, is] = sort(mmm);
bar(log([ide_strength_J(:,is)' ide_strength_J_prior(:,is)']))
set(gca,'xlim',[0 nparam+1])
set(gca,'xticklabel','')
dy = get(gca,'ylim');
for ip=1:nparam,
text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
if iload,
disp('Testing prior mean')
disp_identification(idehess_prior.params, idemodel_prior, idemoments_prior, name);
if advanced, disp('Press ENTER to display advanced diagnostics'), pause, end
plot_identification(idehess_prior.params,idemoments_prior,idehess_prior,idemodel_prior,idelre_prior,advanced,'Prior mean - ',name);
end
legend('relative to param value','relative to prior std','Location','Best')
if flag_score,
title('Identification strength using asymptotic Information matrix (log-scale)')
else
title('Identification strength in the moments (log-scale)')
end
subplot(212)
if SampleSize > 1,
mmm = mean(siJnorm);
else
mmm = (siJnorm);
end
bar(mmm(is))
set(gca,'xlim',[0 nparam+1])
set(gca,'xticklabel','')
dy = get(gca,'ylim');
for ip=1:nparam,
text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
title('Sensitivity in the moments')
if advanced,
figure('Name','Identification in the model'),
subplot(211)
if SampleSize > 1,
mmm = mean(siLREnorm);
else
mmm = (siLREnorm);
end
mmm = [NaN(1, offset), mmm];
bar(mmm(is))
set(gca,'xticklabel','')
for ip=1:nparam,
text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
title('Sensitivity in the LRE model')
subplot(212)
if SampleSize>1,
mmm = mean(siHnorm);
else
mmm = (siHnorm);
end
bar(mmm(is))
set(gca,'xticklabel','')
dy = get(gca,'ylim');
for ip=1:nparam,
text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
title('Sensitivity in the model')
if options_.nograph, close(gcf); end
% identificaton patterns
for j=1:size(cosnJ,2),
pax=NaN(nparam,nparam);
fprintf('\n\n')
disp(['Collinearity patterns with ', int2str(j) ,' parameter(s)'])
fprintf('%-15s [%-*s] %10s\n','Parameter',(15+1)*j,' Expl. params ','cosn')
for i=1:nparam,
namx='';
for in=1:j,
namx=[namx ' ' sprintf('%-15s',name{pars{i,j}(in)})];
fprintf('\n\n')
disp('Testing MC sample')
disp_identification(pdraws, idemodel, idemoments, name);
plot_identification(pdraws,idemoments,idehess_prior,idemodel,idelre,advanced,'MC sample - ',name, IdentifDirectoryName, 1);
if advanced,
jcrit=find(idemoments.ino);
for j=1:length(jcrit),
if ~iload,
[idehess_(j), idemoments_(j), idemodel_(j), idelre_(j), derivatives_info_(j)] = ...
identification_analysis(pdraws(jcrit(j),:),indx,indexo,options_ident,data_info, prior_exist, name_tex,1);
end
fprintf('%-15s [%s] %10.3f\n',name{i},namx,cosnJ(i,j))
pax(i,pars{i,j})=cosnJ(i,j);
end
figure('name',['Collinearity patterns with ', int2str(j) ,' parameter(s)']),
imagesc(pax,[0 1]);
set(gca,'xticklabel','')
set(gca,'yticklabel','')
for ip=1:nparam,
text(ip,(0.5),name{ip},'rotation',90,'HorizontalAlignment','left','interpreter','none')
text(0.5,ip,name{ip},'rotation',0,'HorizontalAlignment','right','interpreter','none')
end
colorbar;
ax=colormap;
ax(1,:)=[0.9 0.9 0.9];
colormap(ax);
if nparam>10,
set(gca,'xtick',(5:5:nparam))
set(gca,'ytick',(5:5:nparam))
end
set(gca,'xgrid','on')
set(gca,'ygrid','on')
saveas(gcf,[IdentifDirectoryName,'/',M_.fname,'_ident_collinearity_', int2str(j)])
eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_collinearity_', int2str(j)]);
eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_collinearity_', int2str(j)]);
if options_.nograph, close(gcf); end
end
disp('')
if flag_score,
[U,S,V]=svd(AHess,0);
if nparam<5,
f1 = figure('name','Identification patterns (Information matrix)');
else
f1 = figure('name','Identification patterns (Information matrix): SMALLEST SV');
f2 = figure('name','Identification patterns (Information matrix): HIGHEST SV');
end
else
[U,S,V]=svd(siJ./normJ(:,ones(nparam,1)),0);
if nparam<5,
f1 = figure('name','Identification patterns (moments)');
else
f1 = figure('name','Identification patterns (moments): SMALLEST SV');
f2 = figure('name','Identification patterns (moments): HIGHEST SV');
tittxt = ['Critical draw n. ',int2str(j),' - '];
fprintf('\n\n')
disp(['Testing ',tittxt, 'Press ENTER']), pause,
plot_identification(pdraws(jcrit(j),:),idemoments_(j),idehess_(j),idemodel_(j),idelre_(j),1,tittxt,name);
% disp_identification(pdraws(jcrit(j),:), idemodel_(j), idemoments_(j), name, advanced);
end
if ~iload,
save([IdentifDirectoryName '/' M_.fname '_identif.mat'], 'idehess_', 'idemoments_','idemodel_', 'idelre_', 'jcrit', '-append');
end
end
for j=1:min(nparam,8),
if j<5,
figure(f1),
jj=j;
else
figure(f2),
jj=j-4;
end
subplot(4,1,jj),
if j<5
bar(abs(V(:,end-j+1))),
Stit = S(end-j+1,end-j+1);
else
bar(abs(V(:,jj))),
Stit = S(jj,jj);
end
set(gca,'xticklabel','')
if j==4 || j==nparam || j==8,
for ip=1:nparam,
text(ip,-0.02,name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
end
title(['Singular value ',num2str(Stit)])
end
figure(f1);
saveas(f1,[IdentifDirectoryName,'/',M_.fname,'_ident_pattern_1'])
eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_pattern_1']);
eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_pattern_1']);
if nparam>4,
figure(f2),
saveas(f2,[IdentifDirectoryName,'/',M_.fname,'_ident_pattern_2'])
eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_pattern_2']);
eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_pattern_2']);
end
if SampleSize>1,
options_.nograph=1;
figure('Name','Condition Number'),
subplot(221)
hist(log10(idemodel.cond))
title('log10 of Condition number in the model')
subplot(222)
hist(log10(idemoments.cond))
title('log10 of Condition number in the moments')
subplot(223)
hist(log10(idelre.cond))
title('log10 of Condition number in the LRE model')
saveas(gcf,[IdentifDirectoryName,'/',M_.fname,'_ident_COND'])
eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_COND']);
eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_COND']);
if options_.nograph, close(gcf); end
ncut=floor(SampleSize/10*9);
[~,is]=sort(idelre.cond);
[proba, dproba] = stab_map_1(pdraws, is(1:ncut), is(ncut+1:end), 'HighestCondNumberLRE', 1, [], IdentifDirectoryName);
[~,is]=sort(idemodel.cond);
[proba, dproba] = stab_map_1(pdraws, is(1:ncut), is(ncut+1:end), 'HighestCondNumberModel', 1, [], IdentifDirectoryName);
[~,is]=sort(idemoments.cond);
[proba, dproba] = stab_map_1(pdraws, is(1:ncut), is(ncut+1:end), 'HighestCondNumberMoments', 1, [], IdentifDirectoryName);
% [proba, dproba] = stab_map_1(idemoments.Mco', is(1:ncut), is(ncut+1:end), 'HighestCondNumberMoments_vs_Mco', 1, [], IdentifDirectoryName);
for j=1:nparam,
% ibeh=find(idemoments.Mco(j,:)<0.9);
% inonbeh=find(idemoments.Mco(j,:)>=0.9);
% if ~isempty(ibeh) && ~isempty(inonbeh)
% [proba, dproba] = stab_map_1(pdraws, ibeh, inonbeh, ['HighestMultiCollinearity_',name{j}], 1, [], IdentifDirectoryName);
% end
[~,is]=sort(idemoments.Mco(j,:));
[proba, dproba] = stab_map_1(pdraws, is(1:ncut), is(ncut+1:end), ['HighestMultiCollinearity_',name{j}], 1, [], IdentifDirectoryName);
end
end
end
if exist('OCTAVE_VERSION')
warning('on'),
else

14
matlab/ident_bruteforce.m
View File

@ -46,11 +46,21 @@ for ll = 1:n,
tmp = find([1:k]~=ii);
tmp2 = nchoosek(tmp,ll);
cosnJ2=zeros(size(tmp2,1),1);
b=[];
for jj = 1:size(tmp2,1)
cosnJ2(jj,1) = cosn([J(:,ii),J(:,tmp2(jj,:))]);
[cosnJ2(jj,1), b(:,jj)] = cosn([J(:,ii),J(:,tmp2(jj,:))]);
end
cosnJ(ii,ll) = max(cosnJ2(:,1));
pars{ii,ll} = tmp2(find(cosnJ2(:,1)==max(cosnJ2(:,1))),:);
if cosnJ(ii,ll)>1.e-8,
if ll>1 && ((cosnJ(ii,ll)-cosnJ(ii,ll-1))<1.e-8),
pars{ii,ll} = [pars{ii,ll-1} NaN];
cosnJ(ii,ll) = cosnJ(ii,ll-1);
else
pars{ii,ll} = tmp2(find(cosnJ2(:,1)==max(cosnJ2(:,1))),:);
end
else
pars{ii,ll} = NaN(1,ll);
end
waitbar(ii/k,h)
end
close(h),

219
matlab/identification_analysis.m Normal file
View File

@ -0,0 +1,219 @@
function [ide_hess, ide_moments, ide_model, ide_lre, derivatives_info, info] = identification_analysis(params,indx,indexo,options_ident,data_info, prior_exist, name_tex, init)
% function [ide_hess, ide_moments, ide_model, ide_lre, derivatives_info, info] = identification_analysis(params,indx,indexo,options_ident,data_info, prior_exist, name_tex, init)
% given the parameter vector params, wraps all identification analyses
%
% INPUTS
% o params [array] parameter values for identification checks
% o indx [array] index of estimated parameters
% o indexo [array] index of estimated shocks
% o options_ident [structure] identification options
% o data_info [structure] data info for Kalmna Filter
% o prior_exist [integer]
% =1 when prior exists and indentification is checked only for estimated params and shocks
% =0 when prior is not defined and indentification is checked for all params and shocks
% o nem_tex [char] list of tex names
% o init [integer] flag for initialization of persistent vars
%
% OUTPUTS
% o ide_hess [structure] identification results using Asymptotic Hessian
% o ide_moments [structure] identification results using theoretical moments
% o ide_model [structure] identification results using reduced form solution
% o ide_lre [structure] identification results using LRE model
% o derivatives_info [structure] info about analytic derivs
% o info output from dynare resolve
%
% SPECIAL REQUIREMENTS
% None
% Copyright (C) 2008-2011 Dynare Team
%
% 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/>.
global oo_ M_ options_ bayestopt_ estim_params_
persistent indH indJJ indLRE
nparam=length(params);
np=length(indx);
offset=nparam-np;
set_all_parameters(params);
nlags = options_ident.ar;
useautocorr = options_ident.useautocorr;
advanced = options_ident.advanced;
replic = options_ident.replic;
periods = options_ident.periods;
max_dim_cova_group = options_ident.max_dim_cova_group;
[I,J]=find(M_.lead_lag_incidence');
ide_hess = struct();
ide_moments = struct();
ide_model = struct();
ide_lre = struct();
derivatives_info = struct();
[A,B,ys,info]=dynare_resolve;
if info(1)==0,
oo0=oo_;
tau=[oo_.dr.ys(oo_.dr.order_var); vec(A); dyn_vech(B*M_.Sigma_e*B')];
yy0=oo_.dr.ys(I);
[residual, g1 ] = feval([M_.fname,'_dynamic'],yy0, ...
oo_.exo_steady_state', M_.params, ...
oo_.dr.ys, 1);
vg1 = [oo_.dr.ys(oo_.dr.order_var); vec(g1)];
[JJ, H, gam, gp, dA, dOm, dYss] = getJJ(A, B, M_,oo0,options_,0,indx,indexo,bayestopt_.mf2,nlags,useautocorr);
derivatives_info.DT=dA;
derivatives_info.DOm=dOm;
derivatives_info.DYss=dYss;
if init,
indJJ = (find(max(abs(JJ'))>1.e-8));
indH = (find(max(abs(H'))>1.e-8));
indLRE = (find(max(abs(gp'))>1.e-8));
end
ide_moments.indJJ=indJJ;
ide_model.indH=indH;
ide_lre.indLRE=indLRE;
TAU(:,1)=tau(indH);
LRE(:,1)=vg1(indLRE);
GAM(:,1)=gam(indJJ);
siJ = (JJ(indJJ,:));
siH = (H(indH,:));
siLRE = (gp(indLRE,:));
normH = max(abs(siH)')';
normJ = max(abs(siJ)')';
normLRE = max(abs(siLRE)')';
ide_moments.siJ=siJ;
ide_model.siH=siH;
ide_lre.siLRE=siLRE;
ide_moments.GAM=GAM;
ide_model.TAU=TAU;
ide_lre.LRE=LRE;
% [ide_checks.idemodel_Mco, ide_checks.idemoments_Mco, ide_checks.idelre_Mco, ...
% ide_checks.idemodel_Pco, ide_checks.idemoments_Pco, ide_checks.idelre_Pco, ...
% ide_checks.idemodel_cond, ide_checks.idemoments_cond, ide_checks.idelre_cond, ...
% ide_checks.idemodel_ee, ide_checks.idemoments_ee, ide_checks.idelre_ee, ...
% ide_checks.idemodel_ind, ide_checks.idemoments_ind, ...
% ide_checks.idemodel_indno, ide_checks.idemoments_indno, ...
% ide_checks.idemodel_ino, ide_checks.idemoments_ino] = ...
% identification_checks(H(indH,:)./normH(:,ones(nparam,1)),JJ(indJJ,:)./normJ(:,ones(nparam,1)), gp(indLRE,:)./normLRE(:,ones(size(gp,2),1)));
[ide_moments.cond, ide_moments.ind0, ide_moments.indno, ide_moments.ino, ide_moments.Mco, ide_moments.Pco, ide_moments.jweak, ide_moments.jweak_pair] = ...
identification_checks(JJ(indJJ,:)./normJ(:,ones(nparam,1)), 0);
[ide_model.cond, ide_model.ind0, ide_model.indno, ide_model.ino, ide_model.Mco, ide_model.Pco, ide_model.jweak, ide_model.jweak_pair] = ...
identification_checks(H(indH,:)./normH(:,ones(nparam,1)), 0);
[ide_lre.cond, ide_lre.ind0, ide_lre.indno, ide_lre.ino, ide_lre.Mco, ide_lre.Pco, ide_lre.jweak, ide_lre.jweak_pair] = ...
identification_checks(gp(indLRE,:)./normLRE(:,ones(size(gp,2),1)), 0);
indok = find(max(ide_moments.indno,[],1)==0);
ide_strength_J=NaN(1,nparam);
ide_strength_J_prior=NaN(1,nparam);
if advanced,
[ide_moments.pars, ide_moments.cosnJ] = ident_bruteforce(JJ(indJJ,:)./normJ(:,ones(nparam,1)),max_dim_cova_group,options_.TeX,name_tex);
end
if init, %~isempty(indok),
normaliz = abs(params);
if prior_exist,
if ~isempty(estim_params_.var_exo),
normaliz1 = estim_params_.var_exo(:,7); % normalize with prior standard deviation
else
normaliz1=[];
end
if ~isempty(estim_params_.param_vals),
normaliz1 = [normaliz1; estim_params_.param_vals(:,7)]'; % normalize with prior standard deviation
end
% normaliz = max([normaliz; normaliz1]);
normaliz1(isinf(normaliz1)) = 1;
else
normaliz1 = NaN(1,nparam);
end
try,
options_.irf = 0;
options_.noprint = 1;
options_.order = 1;
options_.periods = data_info.gend+100;
options_.kalman_algo = 1;
info = stoch_simul(options_.varobs);
datax=oo_.endo_simul(options_.varobs_id,100+1:end);
% datax=data;
derivatives_info.no_DLIK=1;
[fval,cost_flag,ys,trend_coeff,info,DLIK,AHess] = DsgeLikelihood(params',data_info.gend,datax,data_info.data_index,data_info.number_of_observations,data_info.no_more_missing_observations,derivatives_info);
AHess=-AHess;
ide_hess.AHess= AHess;
deltaM = sqrt(diag(AHess));
iflag=any((deltaM.*deltaM)==0);
tildaM = AHess./((deltaM)*(deltaM'));
if iflag || rank(AHess)>rank(tildaM),
[ide_hess.cond, ide_hess.ind0, ide_hess.indno, ide_hess.ino, ide_hess.Mco, ide_hess.Pco] = identification_checks(AHess, 1);
else
[ide_hess.cond, ide_hess.ind0, ide_hess.indno, ide_hess.ino, ide_hess.Mco, ide_hess.Pco] = identification_checks(tildaM, 1);
end
indok = find(max(ide_hess.indno,[],1)==0);
cparam(indok,indok) = inv(AHess(indok,indok));
normaliz(indok) = sqrt(diag(cparam(indok,indok)))';
cmm = siJ*((AHess)\siJ');
rhoM=sqrt(1./diag(inv(tildaM(indok,indok))));
deltaM = deltaM.*params';
flag_score=1;
catch,
replic = max([replic, length(indJJ)*3]);
cmm = simulated_moment_uncertainty(indJJ, periods, replic);
% MIM=siJ(:,indok)'*(cmm\siJ(:,indok));
MIM=siJ'*(cmm\siJ);
ide_hess.AHess= MIM;
deltaM = sqrt(diag(MIM));
iflag=any((deltaM.*deltaM)==0);
tildaM = MIM./((deltaM)*(deltaM'));
if iflag || rank(MIM)>rank(tildaM),
[ide_hess.cond, ide_hess.ind0, ide_hess.indno, ide_hess.ino, ide_hess.Mco, ide_hess.Pco] = identification_checks(MIM, 1);
else
[ide_hess.cond, ide_hess.ind0, ide_hess.indno, ide_hess.ino, ide_hess.Mco, ide_hess.Pco] = identification_checks(tildaM, 1);
end
indok = find(max(ide_hess.indno,[],1)==0);
% rhoM=sqrt(1-1./diag(inv(tildaM)));
% rhoM=(1-1./diag(inv(tildaM)));
rhoM(indok)=sqrt(1./diag(inv(tildaM(indok,indok))));
normaliz(indok) = (sqrt(diag(inv(tildaM(indok,indok))))./deltaM(indok))'; %sqrt(diag(inv(MIM(indok,indok))))';
deltaM = deltaM.*params';
flag_score=0;
end
ide_strength_J(indok) = (1./(normaliz(indok)'./abs(params(indok)')));
ide_strength_J_prior(indok) = (1./(normaliz(indok)'./normaliz1(indok)'));
ide_strength_J(params==0)=ide_strength_J_prior(params==0);
quant = siJ./repmat(sqrt(diag(cmm)),1,nparam);
siJnorm = vnorm(quant).*normaliz1;
% siJnorm = vnorm(siJ(inok,:)).*normaliz;
quant=[];
inok = find((abs(TAU)<1.e-8));
isok = find((abs(TAU)>=1.e-8));
quant(isok,:) = siH(isok,:)./repmat(TAU(isok,1),1,nparam);
quant(inok,:) = siH(inok,:)./repmat(mean(abs(TAU)),length(inok),nparam);
siHnorm = vnorm(quant).*normaliz1;
% siHnorm = vnorm(siH./repmat(TAU,1,nparam)).*normaliz;
quant=[];
inok = find((abs(LRE)<1.e-8));
isok = find((abs(LRE)>=1.e-8));
quant(isok,:) = siLRE(isok,:)./repmat(LRE(isok,1),1,np);
quant(inok,:) = siLRE(inok,:)./repmat(mean(abs(LRE)),length(inok),np);
siLREnorm = vnorm(quant).*normaliz1(offset+1:end);
% siLREnorm = vnorm(siLRE./repmat(LRE,1,nparam-offset)).*normaliz(offset+1:end);
ide_hess.ide_strength_J=ide_strength_J;
ide_hess.ide_strength_J_prior=ide_strength_J_prior;
ide_moments.siJnorm=siJnorm;
ide_model.siHnorm=siHnorm;
ide_lre.siLREnorm=siLREnorm;
ide_hess.flag_score=flag_score;
end,
end

196
matlab/identification_checks.m
View File

@ -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))+1:end) = e1; % non-zero eigenvectors
eH(find(vnorm(H)<eps),1:length(find(vnorm(H)<eps)))=eye(length(find(vnorm(H)<eps)));
condH = cond(H1);
rankH = rank(H);
rankHH = rank(H'*H);
ind2 = find(vnorm(JJ)>=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))+1:end) = ee1; % non-zero eigenvectors
eJ(find(vnorm(JJ)<eps),1:length(find(vnorm(JJ)<eps)))=eye(length(find(vnorm(JJ)<eps)));
condJ = cond(JJ1);
rankJJ = rank(JJ'*JJ);
condJ= cond(JJ1);
rankJ = rank(JJ);
rankJJ = rankJ;
if hess_flag==0,
rankJJ = rank(JJ'*JJ);
end
ind3 = find(vnorm(gp)>=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)<eps))+1:end) = ex1; % non-zero eigenvectors
eGP(find(vnorm(gp)<eps),1:length(find(vnorm(gp)<eps)))=eye(length(find(vnorm(gp)<eps)));
% condJ = cond(JJ1'*JJ1);
condGP = cond(gp1);
ind0 = zeros(1,npar);
ind0(ind1) = 1;
ind01 = zeros(npar,1);
ind02 = zeros(npar,1);
ind01(ind1) = 1;
ind02(ind2) = 1;
% find near linear dependence problems:
McoH = NaN(npar,1);
McoJ = NaN(npar,1);
McoGP = NaN(npar0,1);
for ii = 1:size(H1,2);
McoH(ind1(ii),:) = [cosn([H1(:,ii),H1(:,find([1:1:size(H1,2)]~=ii))])];
end
for ii = 1:size(JJ1,2);
McoJ(ind2(ii),:) = [cosn([JJ1(:,ii),JJ1(:,find([1:1:size(JJ1,2)]~=ii))])];
end
for ii = 1:size(gp1,2);
McoGP(ind3(ii),:) = [cosn([gp1(:,ii),gp1(:,find([1:1:size(gp1,2)]~=ii))])];
end
ixno = 0;
if rankH<npar || rankHH<npar || min(1-McoH)<1.e-10
% - find out which parameters are involved,
% using the vnorm and the svd of H computed before;
% disp('Some parameters are NOT identified in the model: H rank deficient')
% disp(' ')
if length(ind1)<npar,
% parameters with zero column in H
ixno = ixno + 1;
indnoH(ixno,:) = (~ismember([1:npar],ind1));
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))]);
end
e0 = [rankHH+1:length(ind1)];
for j=1:length(e0),
% linearely dependent parameters in H
ixno = ixno + 1;
indnoH(ixno,ind1) = (abs(e1(:,e0(j))) > 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<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,
if length(ind1)<npar,
% parameters with zero column in JJ
ixno = ixno + 1;
indnoJ(ixno,:) = (~ismember([1:npar],ind2));
ixnoJ = ixnoJ + 1;
indnoJ(ixnoJ,:) = (~ismember([1:npar],ind1));
end
ee0 = [rankJJ+1:length(ind2)];
ee0 = [rankJJ+1:length(ind1)];
for j=1:length(ee0),
% linearely dependent parameters in JJ
ixno = ixno + 1;
indnoJ(ixno,ind2) = (abs(ee1(:,ee0(j))) > 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)';

264
matlab/plot_identification.m Normal file
View File

@ -0,0 +1,264 @@
function plot_identification(params,idemoments,idehess,idemodel, idelre, advanced, tittxt, name, IdentifDirectoryName, save_figure)
% function plot_identification(params,idemoments,idehess,idemodel, idelre, advanced, tittxt, name, IdentifDirectoryName, save_figure)
%
% INPUTS
% o params [array] parameter values for identification checks
% o idemoments [structure] identification results for the moments
% o idehess [structure] identification results for the Hessian
% o idemodel [structure] identification results for the reduced form solution
% o idelre [structure] identification results for the LRE model
% o advanced [integer] flag for advanced identification checks
% o tittxt [char] name of the results to plot
% o name [char] list of names
% o IdentifDirectoryName [char] directory name
% o save_figure [integer] flag for saving plots (=1) or not (=0)
%
% OUTPUTS
% None
%
% SPECIAL REQUIREMENTS
% None
% Copyright (C) 2008-2011 Dynare Team
%
% 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/>.
global M_ options_
if nargin<10 || isempty(save_figure),
save_figure=0;
end
[SampleSize, nparam]=size(params);
siJnorm = idemoments.siJnorm;
siHnorm = idemodel.siHnorm;
siLREnorm = idelre.siLREnorm;
% if prior_exist,
% tittxt = 'Prior mean - ';
% else
% tittxt = '';
% end
if SampleSize == 1,
siJ = idemoments.siJ;
normJ = max(abs(siJ)')';
figure('Name',[tittxt, 'Identification using info from observables']),
subplot(211)
mmm = (idehess.ide_strength_J);
[ss, is] = sort(mmm);
bar(log([idehess.ide_strength_J(:,is)' idehess.ide_strength_J_prior(:,is)']))
set(gca,'xlim',[0 nparam+1])
set(gca,'xticklabel','')
dy = get(gca,'ylim');
for ip=1:nparam,
text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
legend('relative to param value','relative to prior std','Location','Best')
if idehess.flag_score,
title('Identification strength in the asymptotic Information matrix (log-scale)')
else
title('Identification strength in the moments (log-scale)')
end
subplot(212)
mmm = (siJnorm)'./max(siJnorm);
if advanced,
mmm1 = (siHnorm)'./max(siHnorm);
mmm=[mmm mmm1];
mmm1 = (siLREnorm)'./max(siLREnorm);
offset=length(siHnorm)-length(siLREnorm);
mmm1 = [NaN(offset,1); mmm1];
mmm=[mmm mmm1];
end
bar(mmm(is,:))
set(gca,'xlim',[0 nparam+1])
set(gca,'xticklabel','')
dy = get(gca,'ylim');
for ip=1:nparam,
text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
if advanced,
legend('Moments','Model','LRE model','Location','Best')
end
title('Sensitivity bars')
if advanced
% identificaton patterns
for j=1:size(idemoments.cosnJ,2),
pax=NaN(nparam,nparam);
fprintf('\n\n')
disp(['Collinearity patterns with ', int2str(j) ,' parameter(s)'])
fprintf('%-15s [%-*s] %10s\n','Parameter',(15+1)*j,' Expl. params ','cosn')
for i=1:nparam,
namx='';
for in=1:j,
dumpindx = idemoments.pars{i,j}(in);
if isnan(dumpindx),
namx=[namx ' ' sprintf('%-15s','--')];
else
namx=[namx ' ' sprintf('%-15s',name{dumpindx})];
pax(i,dumpindx)=idemoments.cosnJ(i,j);
end
end
fprintf('%-15s [%s] %10.3f\n',name{i},namx,idemoments.cosnJ(i,j))
end
figure('name',[tittxt,'Collinearity patterns with ', int2str(j) ,' parameter(s)']),
imagesc(pax,[0 1]);
set(gca,'xticklabel','')
set(gca,'yticklabel','')
for ip=1:nparam,
text(ip,(0.5),name{ip},'rotation',90,'HorizontalAlignment','left','interpreter','none')
text(0.5,ip,name{ip},'rotation',0,'HorizontalAlignment','right','interpreter','none')
end
colorbar;
ax=colormap;
ax(1,:)=[0.9 0.9 0.9];
colormap(ax);
if nparam>10,
set(gca,'xtick',(5:5:nparam))
set(gca,'ytick',(5:5:nparam))
end
set(gca,'xgrid','on')
set(gca,'ygrid','on')
if save_figure
saveas(gcf,[IdentifDirectoryName,'/',M_.fname,'_ident_collinearity_', int2str(j)])
eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_collinearity_', int2str(j)]);
eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_collinearity_', int2str(j)]);
if options_.nograph, close(gcf); end
end
end
disp('')
if idehess.flag_score,
[U,S,V]=svd(idehess.AHess,0);
if nparam<5,
f1 = figure('name',[tittxt,'Identification patterns (Information matrix)']);
else
f1 = figure('name',[tittxt,'Identification patterns (Information matrix): SMALLEST SV']);
f2 = figure('name',[tittxt,'Identification patterns (Information matrix): HIGHEST SV']);
end
else
[U,S,V]=svd(siJ./normJ(:,ones(nparam,1)),0);
if nparam<5,
f1 = figure('name',[tittxt,'Identification patterns (moments)']);
else
f1 = figure('name',[tittxt,'Identification patterns (moments): SMALLEST SV']);
f2 = figure('name',[tittxt,'Identification patterns (moments): HIGHEST SV']);
end
end
for j=1:min(nparam,8),
if j<5,
figure(f1),
jj=j;
else
figure(f2),
jj=j-4;
end
subplot(4,1,jj),
if j<5
bar(abs(V(:,end-j+1))),
Stit = S(end-j+1,end-j+1);
else
bar(abs(V(:,jj))),
Stit = S(jj,jj);
end
set(gca,'xticklabel','')
if j==4 || j==nparam || j==8,
for ip=1:nparam,
text(ip,-0.02,name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
end
title(['Singular value ',num2str(Stit)])
end
if save_figure,
figure(f1);
saveas(f1,[IdentifDirectoryName,'/',M_.fname,'_ident_pattern_1'])
eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_pattern_1']);
eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_pattern_1']);
if nparam>4,
figure(f2),
saveas(f2,[IdentifDirectoryName,'/',M_.fname,'_ident_pattern_2'])
eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_pattern_2']);
eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_pattern_2']);
end
end
end
else
figure('Name',['MC sensitivities']),
mmm = (idehess.ide_strength_J);
[ss, is] = sort(mmm);
mmm = mean(siJnorm)';
mmm = mmm./max(mmm);
if advanced,
mmm1 = mean(siHnorm)';
mmm=[mmm mmm1./max(mmm1)];
mmm1 = mean(siLREnorm)';
offset=size(siHnorm,2)-size(siLREnorm,2);
mmm1 = [NaN(offset,1); mmm1./max(mmm1)];
mmm=[mmm mmm1];
end
bar(mmm(is,:))
set(gca,'xlim',[0 nparam+1])
set(gca,'xticklabel','')
dy = get(gca,'ylim');
for ip=1:nparam,
text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
if advanced,
legend('Moments','Model','LRE model','Location','Best')
end
title('MC mean of sensitivity measures')
if advanced,
options_.nograph=1;
figure('Name','MC Condition Number'),
subplot(221)
hist(log10(idemodel.cond))
title('log10 of Condition number in the model')
subplot(222)
hist(log10(idemoments.cond))
title('log10 of Condition number in the moments')
subplot(223)
hist(log10(idelre.cond))
title('log10 of Condition number in the LRE model')
saveas(gcf,[IdentifDirectoryName,'/',M_.fname,'_ident_COND'])
eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_COND']);
eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_COND']);
if options_.nograph, close(gcf); end
ncut=floor(SampleSize/10*9);
[~,is]=sort(idelre.cond);
[proba, dproba] = stab_map_1(params, is(1:ncut), is(ncut+1:end), 'MC_HighestCondNumberLRE', 1, [], IdentifDirectoryName, 0.1);
[~,is]=sort(idemodel.cond);
[proba, dproba] = stab_map_1(params, is(1:ncut), is(ncut+1:end), 'MC_HighestCondNumberModel', 1, [], IdentifDirectoryName, 0.1);
[~,is]=sort(idemoments.cond);
[proba, dproba] = stab_map_1(params, is(1:ncut), is(ncut+1:end), 'MC_HighestCondNumberMoments', 1, [], IdentifDirectoryName, 0.1);
% [proba, dproba] = stab_map_1(idemoments.Mco', is(1:ncut), is(ncut+1:end), 'HighestCondNumberMoments_vs_Mco', 1, [], IdentifDirectoryName);
for j=1:nparam,
% ibeh=find(idemoments.Mco(j,:)<0.9);
% inonbeh=find(idemoments.Mco(j,:)>=0.9);
% if ~isempty(ibeh) && ~isempty(inonbeh)
% [proba, dproba] = stab_map_1(params, ibeh, inonbeh, ['HighestMultiCollinearity_',name{j}], 1, [], IdentifDirectoryName);
% end
[~,is]=sort(idemoments.Mco(:,j));
[proba, dproba] = stab_map_1(params, is(1:ncut), is(ncut+1:end), ['MC_HighestMultiCollinearity_',name{j}], 1, [], IdentifDirectoryName, 0.15);
end
end
end
% disp_identification(params, idemodel, idemoments, name)