Modifications following meetings in Paris (not yet completed 100%):

1) new default option for prior_mc = 1;
2) new option: advanced (default=0). If advanced=0, simple messages and plots; advanced=1 shows all details about weak identification diagnostics;
3) new option: periods (default=300). Length of stochastic simulations to generate simulated moments for 'GMM type' analysis of strength of identification;
4) new option: replic (default = 100). Number of replicated simulated moments for 'GMM type' analysis of strength of identification;
5) new feature: allow local identification analysis WITHOUT defining priors, using parameter values as defined in the .mod file definition.
6) simplified graphical and tabular output.
time-shift
Marco Ratto 2010-09-29 16:12:48 +02:00
parent 0a637be94c
commit 70400e5e84
1 changed files with 475 additions and 144 deletions

View File

@ -30,14 +30,27 @@ end
options_ident = set_default_option(options_ident,'load_ident_files',0);
options_ident = set_default_option(options_ident,'useautocorr',0);
options_ident = set_default_option(options_ident,'ar',3);
options_ident = set_default_option(options_ident,'prior_mc',2000);
options_ident = set_default_option(options_ident,'prior_mc',1);
options_ident = set_default_option(options_ident,'prior_range',0);
options_ident = set_default_option(options_ident,'periods',300);
options_ident = set_default_option(options_ident,'replic',100);
options_ident = set_default_option(options_ident,'advanced',0);
if nargin==2,
options_ident.prior_mc=size(pdraws0,1);
end
if isempty(estim_params_),
options_ident.prior_mc=1;
options_ident.prior_range=0;
prior_exist=0;
else
prior_exist=1;
end
iload = options_ident.load_ident_files;
advanced = options_ident.advanced;
nlags = options_ident.ar;
periods = options_ident.periods;
replic = options_ident.replic;
useautocorr = options_ident.useautocorr;
options_.ar=nlags;
options_.prior_mc = options_ident.prior_mc;
@ -68,33 +81,60 @@ if ~(exist('sylvester3mr','file')==2),
end
IdentifDirectoryName = CheckPath('identification');
if prior_exist,
indx = estim_params_.param_vals(:,1);
indexo=[];
if ~isempty(estim_params_.var_exo)
indexo = estim_params_.var_exo(:,1);
indx = estim_params_.param_vals(:,1);
indexo=[];
if ~isempty(estim_params_.var_exo)
indexo = estim_params_.var_exo(:,1);
end
nparam = length(bayestopt_.name);
np = estim_params_.np;
name = bayestopt_.name;
offset = estim_params_.nvx;
offset = offset + estim_params_.nvn;
offset = offset + estim_params_.ncx;
offset = offset + estim_params_.ncn;
else
indx = [1:M_.param_nbr];
indexo = [1:M_.exo_nbr];
offset = M_.exo_nbr;
np = M_.param_nbr;
nparam = np+offset;
name = [cellstr(M_.exo_names); cellstr(M_.param_names)];
end
nparam = length(bayestopt_.name);
MaxNumberOfBytes=options_.MaxNumberOfBytes;
disp(' ')
disp(['==== Identification analysis ====' ]),
disp(' ')
if iload <=0,
iteration = 0;
burnin_iteration = 0;
if SampleSize==1,
BurninSampleSize=0;
else
BurninSampleSize=50;
end
loop_indx = 0;
file_index = 0;
run_index = 0;
if SampleSize > 1,
h = waitbar(0,'Monte Carlo identification checks ...');
end
[I,J]=find(M_.lead_lag_incidence');
while iteration < SampleSize,
loop_indx = loop_indx+1;
if prior_exist,
if nargin==2,
if burnin_iteration>=50,
if burnin_iteration>=BurninSampleSize,
params = pdraws0(iteration+1,:);
else
params = pdraws0(burnin_iteration+1,:);
@ -103,6 +143,9 @@ if iload <=0,
params = prior_draw();
end
set_all_parameters(params);
else
params = [sqrt(diag(M_.Sigma_e))', M_.params'];
end
[A,B,ys,info]=dynare_resolve;
@ -117,7 +160,7 @@ if iload <=0,
yy0=oo_.dr.ys(I);
[residual, g1 ] = feval([M_.fname,'_dynamic'],yy0, oo_.exo_steady_state', M_.params,1);
if burnin_iteration<50,
if burnin_iteration<BurninSampleSize,
burnin_iteration = burnin_iteration + 1;
pdraws(burnin_iteration,:) = params;
TAU(:,burnin_iteration)=tau;
@ -143,11 +186,11 @@ if iload <=0,
dum = [dum; vec(gam{j+1})];
end
GAM(:,burnin_iteration)=[oo_.dr.ys(bayestopt_.mfys); dum];
warning warning_old_state;
% warning warning_old_state;
else
iteration = iteration + 1;
run_index = run_index + 1;
if iteration==1,
if iteration==1 & BurninSampleSize,
indJJ = (find(std(GAM')>1.e-8));
indH = (find(std(TAU')>1.e-8));
indLRE = (find(std(LRE')>1.e-8));
@ -163,10 +206,23 @@ if iload <=0,
end
if iteration,
[JJ, H, gam, gp] = getJJ(A, B, M_,oo0,options_,0,indx,indexo,bayestopt_.mf2,nlags,useautocorr);
if BurninSampleSize == 0,
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);
[JJ, H, gam, gp] = getJJ(A, B, M_,oo0,options_,0,indx,indexo,bayestopt_.mf2,nlags,useautocorr);
GAM(:,iteration)=gam(indJJ);
stoLRE(:,:,run_index) = gp(indLRE,:);
stoH(:,:,run_index) = H(indH,:);
@ -177,12 +233,16 @@ if iload <=0,
% use prior uncertainty
siJ = (JJ(indJJ,:));
siH = (H(indH,:));
siLRE = (gp(indLRE,:));
% siJ = abs(JJ(indJJ,:).*(ones(length(indJJ),1)*bayestopt_.p2'));
% siH = abs(H(indH,:).*(ones(length(indH),1)*bayestopt_.p2'));
% siJ = abs(JJ(indJJ,:).*(1./mGAM'*bayestopt_.p2'));
% siH = abs(H(indH,:).*(1./mTAU'*bayestopt_.p2'));
siJnorm(iteration,:) = vnorm(siJ./repmat(GAM(:,iteration),1,nparam)).*params;
siHnorm(iteration,:) = vnorm(siH./repmat(TAU(:,iteration),1,nparam)).*params;
siLREnorm(iteration,:) = vnorm(siLRE./repmat(LRE(:,iteration),1,nparam-offset)).*params(offset+1:end);
if iteration ==1,
siJmean = abs(siJ)./SampleSize;
siHmean = abs(siH)./SampleSize;
@ -210,8 +270,36 @@ if iload <=0,
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.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)), bayestopt_);
% identification_checks(H(indH,:),JJ(indJJ,:), gp(indLRE,:), bayestopt_);
indok = find(max(idemoments.indno{iteration},[],1)==0);
ide_strength_J(iteration,:)=NaN(1,nparam);
if iteration ==1,
cmm = simulated_moment_uncertainty(indJJ, periods, replic);
end,
% Jinv=(siJ(:,indok)'*siJ(:,indok))\siJ(:,indok)';
% MIM=inv(Jinv*cmm*Jinv');
MIM=siJ(:,indok)'*inv(cmm)*siJ(:,indok);
deltaM = sqrt(diag(MIM));
tildaM = MIM./((deltaM)*(deltaM'));
rhoM=sqrt(1-1./diag(inv(tildaM)));
deltaM = deltaM.*params(indok)';
ide_strength_J(iteration,indok) = (1./[sqrt(diag(inv(MIM)))./params(indok)']);
% indok = find(max(idemodel.indno{iteration},[],1)==0);
% ide_strength_H(iteration,:)=zeros(1,nparam);
% mim=inv(siH(:,indok)'*siH(:,indok))*siH(:,indok)';
% % mim=mim*diag(GAM(:,iteration))*mim';
% % MIM=inv(mim);
% mim=mim.*repmat(TAU(:,iteration),1,length(indok))';
% MIM=inv(mim*mim');
% deltaM = sqrt(diag(MIM));
% tildaM = MIM./((deltaM)*(deltaM'));
% rhoM=sqrt(1-1./diag(inv(tildaM)));
% deltaM = deltaM.*params(indok)';
% ide_strength_H(iteration,indok) = (1./[sqrt(diag(inv(MIM)))./params(indok)']);
if run_index==MAX_tau | iteration==SampleSize,
file_index = file_index + 1;
if run_index<MAX_tau,
@ -224,13 +312,17 @@ if iload <=0,
end
if SampleSize > 1,
waitbar(iteration/SampleSize,h,['MC Identification checks ',int2str(iteration),'/',int2str(SampleSize)])
end
end
end
end
if SampleSize > 1,
close(h)
end
save([IdentifDirectoryName '/' M_.fname '_identif'], 'pdraws', 'idemodel', 'idemoments', 'idelre', 'indJJ', 'indH', 'indLRE', ...
@ -263,7 +355,8 @@ else
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.indno{iteration}, idemoments.indno{iteration}, ...
idemodel.ino(iteration), idemoments.ino(iteration)] = ...
identification_checks(stoH(:,:,index)./normH(:,ones(nparam,1)), ...
stoJJ(:,:,index)./normJ(:,ones(nparam,1)), ...
stoLRE(:,:,index)./normLRE(:,ones(size(stoLRE,2),1)), bayestopt_);
@ -273,13 +366,28 @@ else
close(h);
save([IdentifDirectoryName '/' M_.fname '_identif'], 'idemodel', 'idemoments', 'idelre', '-append')
end
iteration = 0;
h = waitbar(0,'Monte Carlo identification checks ...');
for file_index=1:length(identFiles)
load([IdentifDirectoryName '/' M_.fname '_identif_' int2str(file_index)], 'stoH', 'stoJJ', 'stoLRE')
for index=1:size(stoH,3),
iteration = iteration+1;
fobj(iteration)=sum(((GAM(:,iteration)-GAM(:,1))).^2);
fobjH(iteration)=sum(((TAU(:,iteration)-TAU(:,1))).^2);
fobjR(iteration)=sum(((GAM(:,iteration)./GAM(:,1)-1)).^2);
fobjHR(iteration)=sum(((TAU(:,iteration)./TAU(:,1)-1)).^2);
FOC(iteration,:) = (GAM(:,iteration)-GAM(:,1))'*stoJJ(:,:,index);
FOCH(iteration,:) = (TAU(:,iteration)-TAU(:,1))'*stoH(:,:,index);
FOCR(iteration,:) = ((GAM(:,iteration)./GAM(:,1)-1)./GAM(:,1))'*stoJJ(:,:,index);
FOCHR(iteration,:) = ((TAU(:,iteration)./TAU(:,1)-1)./TAU(:,1))'*stoH(:,:,index);
waitbar(iteration/SampleSize,h,['MC Identification checks ',int2str(iteration),'/',int2str(SampleSize)])
end
end
close(h);
end
offset = estim_params_.nvx;
offset = offset + estim_params_.nvn;
offset = offset + estim_params_.ncx;
offset = offset + estim_params_.ncn;
if SampleSize>1,
siJmean = siJmean.*(ones(length(indJJ),1)*std(pdraws));
siHmean = siHmean.*(ones(length(indH),1)*std(pdraws));
siLREmean = siLREmean.*(ones(length(indLRE),1)*std(pdraws(:, offset+1:end )));
@ -290,15 +398,16 @@ derLREmean = derLREmean.*(ones(length(indLRE),1)*std(pdraws(:, offset+1:end )));
derHmean = abs(derHmean./(max(siHmean')'*ones(1,size(pdraws,2))));
derJmean = abs(derJmean./(max(siJmean')'*ones(1,size(pdraws,2))));
derLREmean = abs(derLREmean./(max(siLREmean')'*ones(1,estim_params_.np)));
derLREmean = abs(derLREmean./(max(siLREmean')'*ones(1,np)));
siHmean = siHmean./(max(siHmean')'*ones(1,size(pdraws,2)));
siJmean = siJmean./(max(siJmean')'*ones(1,size(pdraws,2)));
siLREmean = siLREmean./(max(siLREmean')'*ones(1,estim_params_.np));
siLREmean = siLREmean./(max(siLREmean')'*ones(1,np));
tstJmean = derJmean*0;
tstHmean = derHmean*0;
tstLREmean = derLREmean*0;
if exist('OCTAVE_VERSION')
warning('off', 'Octave:divide-by-zero')
else
@ -315,6 +424,7 @@ if j>offset
tstLREmean(indd,j-offset) = abs(derLREmean(indd,j-offset))./siLREmean(indd,j-offset);
end
end
end
if nargout>3 & iload,
@ -382,7 +492,7 @@ end
% siHPCA2 = siHPCA2./max(siHPCA2);
disp_identification(pdraws, idemodel, idemoments)
disp_identification(pdraws, idemodel, idemoments, name, advanced)
% figure,
% % myboxplot(siPCA(1:(max(find(cumsum(latent)./length(indJJ)<0.99))+1),:))
@ -392,7 +502,7 @@ disp_identification(pdraws, idemodel, idemoments)
% set(gca,'xticklabel','')
% set(gca,'xlim',[0.5 nparam+0.5])
% for ip=1:nparam,
% text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
% text(ip,-0.02,name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
% end
% title('Sensitivity in TAU''s PCA')
%
@ -403,7 +513,7 @@ disp_identification(pdraws, idemodel, idemoments)
% set(gca,'xticklabel','')
% set(gca,'xlim',[0.5 nparam+0.5])
% for ip=1:nparam,
% text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
% text(ip,-0.02,name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
% end
% title('Sensitivity in standardized TAU''s PCA')
%
@ -415,7 +525,7 @@ disp_identification(pdraws, idemodel, idemoments)
% set(gca,'xticklabel','')
% set(gca,'xlim',[0.5 nparam+0.5])
% for ip=1:nparam,
% text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
% text(ip,-0.02,name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
% end
% title('Sensitivity in moments'' PCA')
%
@ -426,30 +536,50 @@ disp_identification(pdraws, idemodel, idemoments)
% set(gca,'xticklabel','')
% set(gca,'xlim',[0.5 nparam+0.5])
% for ip=1:nparam,
% text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
% text(ip,-0.02,name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
% end
% title('Sensitivity in standardized moments'' PCA')
if advanced,
figure('Name','Identification LRE model form'),
subplot(211)
mmm = mean(siLREmean);
if SampleSize > 1,
mmm = mean(siLREnorm);
else
mmm = (siLREnorm);
end
[ss, is] = sort(mmm);
myboxplot(siLREmean(:,is))
set(gca,'ylim',[0 1.05])
if SampleSize ==1,
bar(siLREnorm(:,is))
else
myboxplot(log(siLREnorm(:,is)))
end
% mmm = mean(siLREmean);
% [ss, is] = sort(mmm);
% myboxplot(siLREmean(:,is))
% set(gca,'ylim',[0 1.05])
set(gca,'xticklabel','')
for ip=1:estim_params_.np,
text(ip,-0.02,deblank(M_.param_names(estim_params_.param_vals(is(ip),1),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
for ip=1:np,
text(ip,-0.02,deblank(M_.param_names(indx(is(ip)),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
title('Sensitivity in the LRE model')
subplot(212)
if SampleSize>1,
mmm = mean(-idelre.Mco');
else
mmm = (-idelre.Mco');
end
[ss, is] = sort(mmm);
myboxplot(idelre.Mco(is,:)')
if SampleSize ==1,
bar(idelre.Mco(is,:)')
else
myboxplot(idelre.Mco(is,:)')
end
set(gca,'ylim',[0 1])
set(gca,'xticklabel','')
for ip=1:estim_params_.np,
text(ip,-0.02,deblank(M_.param_names(estim_params_.param_vals(is(ip),1),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
for ip=1:np,
text(ip,-0.02,deblank(M_.param_names(indx(is(ip)),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
title('Multicollinearity in the LRE model')
saveas(gcf,[IdentifDirectoryName,'/',M_.fname,'_ident_LRE'])
@ -458,60 +588,176 @@ eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_LRE']);
if options_.nograph, close(gcf); end
figure('Name','Identification in the model'),
% subplot(311)
%
% if SampleSize>1,
% mmm = mean(ide_strength_H);
% else
% mmm = (ide_strength_H);
% end
% [ss, is] = sort(mmm);
% if SampleSize>1,
% myboxplot(ide_strength_H(:,is))
% else
% bar(ide_strength_H(:,is))
% end
% % set(gca,'ylim',[0 1.05])
% set(gca,'xticklabel','')
% dy = get(gca,'ylim');
% % dy=dy(2)-dy(1);
% for ip=1:nparam,
% text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
% end
% title('Identification strength in the model')
subplot(211)
mmm = mean(siHmean);
% mmm = mean(siHmean);
% [ss, is] = sort(mmm);
% myboxplot(siHmean(:,is))
if SampleSize>1,
mmm = mean(siHnorm);
else
mmm = (siHnorm);
end
[ss, is] = sort(mmm);
myboxplot(siHmean(:,is))
set(gca,'ylim',[0 1.05])
if SampleSize>1,
myboxplot(log(siHnorm(:,is)))
else
bar(siHnorm(:,is))
end
% set(gca,'ylim',[0 1.05])
set(gca,'xticklabel','')
dy = get(gca,'ylim');
% dy=dy(2)-dy(1);
for ip=1:nparam,
text(ip,-0.02,bayestopt_.name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
title('Sensitivity in the model')
subplot(212)
if SampleSize>1,
mmm = mean(-idemodel.Mco');
[ss, is] = sort(mmm);
else
mmm = (-idemodel.Mco');
end
% [ss, is] = sort(mmm);
if SampleSize>1,
myboxplot(idemodel.Mco(is,:)')
else
bar(idemodel.Mco(is,:)')
end
set(gca,'ylim',[0 1])
set(gca,'xticklabel','')
for ip=1:nparam,
text(ip,-0.02,bayestopt_.name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
text(ip,-0.02,name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
title('Multicollinearity in the model')
saveas(gcf,[IdentifDirectoryName,'/',M_.fname,'_ident_model'])
eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_model']);
eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_model']);
if options_.nograph, close(gcf); end
end
figure('Name','Identification in the moments'),
subplot(211)
mmm = mean(siJmean);
% mmm = mean(siHmean);
% [ss, is] = sort(mmm);
% myboxplot(siHmean(:,is))
if SampleSize>1,
mmm = mean(ide_strength_J);
else
mmm = (ide_strength_J);
end
[ss, is] = sort(mmm);
myboxplot(siJmean(:,is))
set(gca,'ylim',[0 1.05])
if SampleSize>1,
myboxplot(log(ide_strength_J(:,is)))
else
bar(ide_strength_J(:,is))
end
% set(gca,'ylim',[0 1.05])
set(gca,'xticklabel','')
dy = get(gca,'ylim');
% dy=dy(2)-dy(1);
for ip=1:nparam,
text(ip,-0.02,bayestopt_.name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
title('Identification strength in the moments')
subplot(212)
% mmm = mean(siJmean);
% [ss, is] = sort(mmm);
% myboxplot(siJmean(:,is))
if SampleSize > 1,
mmm = mean(siJnorm);
else
mmm = (siJnorm);
end
% [ss, is] = sort(mmm);
if SampleSize > 1,
myboxplot(log(siJnorm(:,is)))
else
bar(siJnorm(:,is))
end
% set(gca,'ylim',[0 1.05])
set(gca,'xticklabel','')
dy = get(gca,'ylim');
% dy=dy(2)-dy(1);
for ip=1:nparam,
text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
title('Sensitivity in the moments')
subplot(212)
mmm = mean(-idemoments.Mco');
[ss, is] = sort(mmm);
myboxplot(idemoments.Mco(is,:)')
set(gca,'ylim',[0 1])
set(gca,'xticklabel','')
for ip=1:nparam,
text(ip,-0.02,bayestopt_.name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
% subplot(313)
% if SampleSize>1,
% mmm = mean(-idemoments.Mco');
% else
% mmm = (-idemoments.Mco');
% end
% % [ss, is] = sort(mmm);
% if SampleSize>1,
% myboxplot(idemoments.Mco(is,:)')
% else
% bar(idemoments.Mco(is,:)')
% end
% set(gca,'ylim',[0 1])
% set(gca,'xticklabel','')
% for ip=1:nparam,
% text(ip,-0.02,name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
% end
% title('Multicollinearity in the moments')
% saveas(gcf,[IdentifDirectoryName,'/',M_.fname,'_ident_moments'])
% eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_moments']);
% eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_moments']);
% if options_.nograph, close(gcf); end
if SampleSize==1,
% identificaton patterns
[U,S,V]=svd(siJ./normJ(:,ones(nparam,1)),0);
figure('name','Identification patterns (moments)'),
for j=1:min(nparam,8),
subplot(2,4,j),
if j<5
bar(abs(V(:,end-j+1))),
Stit = S(end-j+1,end-j+1);
if j==1, ylabel('SMALLEST singular values'), end,
else
bar(abs(V(:,j))),
Stit = S(j,j);
if j==5, ylabel('LARGEST singular values'), end,
end
set(gca,'xticklabel','')
for ip=1:nparam,
text(ip,-0.02,name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
title(['Singular value ',num2str(Stit)])
end
saveas(gcf,[IdentifDirectoryName,'/',M_.fname,'_ident_pattern'])
eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_pattern']);
eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_pattern']);
end
title('Multicollinearity in the moments')
saveas(gcf,[IdentifDirectoryName,'/',M_.fname,'_ident_moments'])
eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_moments']);
eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_moments']);
if options_.nograph, close(gcf); end
if advanced,
if SampleSize>1
figure('Name','Condition Number'),
subplot(221)
hist(log10(idemodel.cond))
@ -526,98 +772,178 @@ 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
ifig=0;
nbox = min(estim_params_.np-1,12);
for j=1:estim_params_.np,
if mod(j,12)==1,
ifig = ifig+1;
figure('name','Partial correlations in the LRE model'),
iplo=0;
end
iplo=iplo+1;
mmm = mean(squeeze(idelre.Pco(:,j,:))');
[sss, immm] = sort(-mmm);
subplot(3,4,iplo),
if nbox==1,
myboxplot(squeeze(idelre.Pco(immm(2:nbox+1),j,:))),
else
myboxplot(squeeze(idelre.Pco(immm(2:nbox+1),j,:))'),
end
set(gca,'ylim',[0 1])
set(gca,'xticklabel','')
for ip=1:nbox, %estim_params_.np,
text(ip,-0.02,deblank(M_.param_names(estim_params_.param_vals(immm(ip+1),1),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
title(deblank(M_.param_names(estim_params_.param_vals(j,1),:))),
if j==estim_params_.np | mod(j,12)==0
saveas(gcf,[IdentifDirectoryName,'/',M_.fname,'_ident_PCORR_LRE',int2str(ifig)])
eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_PCORR_LRE',int2str(ifig)]);
eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_PCORR_LRE',int2str(ifig)]);
if options_.nograph, close(gcf); end
end
end
ifig=0;
nbox = min(nparam-1,12);
for j=1:nparam,
if mod(j,12)==1,
ifig = ifig+1;
figure('name','Partial correlations in the model'),
iplo=0;
end
iplo=iplo+1;
mmm = mean(squeeze(idemodel.Pco(:,j,:))');
[sss, immm] = sort(-mmm);
subplot(3,4,iplo),
if nbox==1,
myboxplot(squeeze(idemodel.Pco(immm(2:nbox+1),j,:))),
else
myboxplot(squeeze(idemodel.Pco(immm(2:nbox+1),j,:))'),
end
set(gca,'ylim',[0 1])
set(gca,'xticklabel','')
for ip=1:nbox, %estim_params_.np,
text(ip,-0.02,bayestopt_.name{immm(ip+1)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
title(bayestopt_.name{j}),
if j==nparam | mod(j,12)==0
saveas(gcf,[IdentifDirectoryName,'/',M_.fname,'_ident_PCORR_model',int2str(ifig)])
eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_PCORR_model',int2str(ifig)]);
eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_PCORR_model',int2str(ifig)]);
if options_.nograph, close(gcf); end
end
end
ifig=0;
nbox = min(nparam-1,12);
for j=1:nparam,
if mod(j,12)==1,
ifig = ifig+1;
figure('name','Partial correlations in the 1st and 2nd moments'),
iplo=0;
end
iplo=iplo+1;
mmm = mean(squeeze(idemoments.Pco(:,j,:))');
[sss, immm] = sort(-mmm);
subplot(3,4,iplo),
if nbox==1,
myboxplot(squeeze(idemoments.Pco(immm(2:nbox+1),j,:))),
pco = NaN(np,np);
for j=1:np,
if SampleSize>1
pco(j+1:end,j) = mean(abs(squeeze(idelre.Pco(j+1:end,j,:))'));
else
myboxplot(squeeze(idemoments.Pco(immm(2:nbox+1),j,:))'),
end
set(gca,'ylim',[0 1])
set(gca,'xticklabel','')
for ip=1:nbox, %estim_params_.np,
text(ip,-0.02,bayestopt_.name{immm(ip+1)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
title(bayestopt_.name{j}),
if j==nparam | mod(j,12)==0
saveas(gcf,[IdentifDirectoryName,'/',M_.fname,'_ident_PCORR_moments',int2str(ifig)])
eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_PCORR_moments',int2str(ifig)]);
eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_PCORR_moments',int2str(ifig)]);
if options_.nograph, close(gcf); end
pco(j+1:end,j) = abs(idelre.Pco(j+1:end,j))';
end
end
figure('name','Pairwise correlations in the LRE model'),
imagesc(pco',[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;
saveas(gcf,[IdentifDirectoryName,'/',M_.fname,'_ident_PCORR_LRE'])
eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_PCORR_LRE']);
eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_PCORR_LRE']);
if options_.nograph, close(gcf); end
pco = NaN(nparam,nparam);
for j=1:nparam,
if SampleSize>1
pco(j+1:end,j) = mean(abs(squeeze(idemodel.Pco(j+1:end,j,:))'));
else
pco(j+1:end,j) = abs(idemodel.Pco(j+1:end,j))';
end
end
figure('name','Pairwise correlations in the model'),
imagesc(pco',[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;
saveas(gcf,[IdentifDirectoryName,'/',M_.fname,'_ident_PCORR_model'])
eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_PCORR_model']);
eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_PCORR_model']);
if options_.nograph, close(gcf); end
for j=1:nparam,
if SampleSize>1
pco(j+1:end,j) = mean(abs(squeeze(idemoments.Pco(j+1:end,j,:))'));
else
pco(j+1:end,j) = abs(idemoments.Pco(j+1:end,j))';
end
end
figure('name','Pairwise correlations in the moments'),
imagesc(pco',[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;
saveas(gcf,[IdentifDirectoryName,'/',M_.fname,'_ident_PCORR_moments'])
eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_PCORR_moments']);
eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_PCORR_moments']);
if options_.nograph, close(gcf); end
end
% ifig=0;
% nbox = min(np-1,12);
% for j=1:np,
% if mod(j,12)==1,
% ifig = ifig+1;
% figure('name','Pairwise correlations in the LRE model'),
% iplo=0;
% end
% iplo=iplo+1;
% if SampleSize>1
% mmm = mean(abs(squeeze(idelre.Pco(:,j,:))'));
% else
% mmm = abs(idelre.Pco(:,j))';
% end
% [sss, immm] = sort(-mmm);
% subplot(3,4,iplo),
% if nbox==1,
% myboxplot(squeeze(idelre.Pco(immm(2:nbox+1),j,:))),
% else
% myboxplot(squeeze(idelre.Pco(immm(2:nbox+1),j,:))'),
% end
% set(gca,'ylim',[-1 1],'ygrid','on')
% set(gca,'xticklabel','')
% for ip=1:nbox, %np,
% text(ip,-1.02,deblank(M_.param_names(indx(immm(ip+1)),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
% end
% title(deblank(M_.param_names(indx(j),:))),
% if j==np | mod(j,12)==0
% saveas(gcf,[IdentifDirectoryName,'/',M_.fname,'_ident_PCORR_LRE',int2str(ifig)])
% eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_PCORR_LRE',int2str(ifig)]);
% eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_PCORR_LRE',int2str(ifig)]);
% if options_.nograph, close(gcf); end
% end
% end
%
% ifig=0;
% nbox = min(nparam-1,12);
% for j=1:nparam,
% if mod(j,12)==1,
% ifig = ifig+1;
% figure('name','Pairwise correlations in the model'),
% iplo=0;
% end
% iplo=iplo+1;
% if SampleSize>1,
% mmm = mean(abs(squeeze(idemodel.Pco(:,j,:))'));
% else
% mmm = abs(idemodel.Pco(:,j)');
% end
% [sss, immm] = sort(-mmm);
% subplot(3,4,iplo),
% if nbox==1,
% myboxplot(squeeze(idemodel.Pco(immm(2:nbox+1),j,:))),
% else
% myboxplot(squeeze(idemodel.Pco(immm(2:nbox+1),j,:))'),
% end
% set(gca,'ylim',[-1 1],'ygrid','on')
% set(gca,'xticklabel','')
% for ip=1:nbox, %np,
% text(ip,-1.02,name{immm(ip+1)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
% end
% title(name{j}),
% if j==nparam | mod(j,12)==0
% saveas(gcf,[IdentifDirectoryName,'/',M_.fname,'_ident_PCORR_model',int2str(ifig)])
% eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_PCORR_model',int2str(ifig)]);
% eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_PCORR_model',int2str(ifig)]);
% if options_.nograph, close(gcf); end
% end
% end
%
% ifig=0;
% nbox = min(nparam-1,12);
% for j=1:nparam,
% if mod(j,12)==1,
% ifig = ifig+1;
% figure('name','Pairwise correlations in the 1st and 2nd moments'),
% iplo=0;
% end
% iplo=iplo+1;
% if SampleSize>1
% mmm = mean(abs(squeeze(idemoments.Pco(:,j,:))'));
% else
% mmm = abs(idemoments.Pco(:,j)');
% end
% [sss, immm] = sort(-mmm);
% subplot(3,4,iplo),
% if nbox==1,
% myboxplot(squeeze(idemoments.Pco(immm(2:nbox+1),j,:))),
% else
% myboxplot(squeeze(idemoments.Pco(immm(2:nbox+1),j,:))'),
% end
% set(gca,'ylim',[-1 1],'ygrid','on')
% set(gca,'xticklabel','')
% for ip=1:nbox, %np,
% text(ip,-1.02,name{immm(ip+1)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
% end
% title(name{j}),
% if j==nparam | mod(j,12)==0
% saveas(gcf,[IdentifDirectoryName,'/',M_.fname,'_ident_PCORR_moments',int2str(ifig)])
% eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_PCORR_moments',int2str(ifig)]);
% eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_PCORR_moments',int2str(ifig)]);
% if options_.nograph, close(gcf); end
% end
% end
if exist('OCTAVE_VERSION')
@ -625,3 +951,8 @@ if exist('OCTAVE_VERSION')
else
warning on MATLAB:dividebyzero
end
disp(' ')
disp(['==== Identification analysis completed ====' ]),
disp(' ')
disp(' ')