diff --git a/matlab/dynare_identification.m b/matlab/dynare_identification.m index 263d347e8..ea0047abf 100644 --- a/matlab/dynare_identification.m +++ b/matlab/dynare_identification.m @@ -1,4 +1,21 @@ function [pdraws, TAU, GAM, LRE, gp, H, JJ] = dynare_identification(options_ident, pdraws0) +%function [pdraws, TAU, GAM, LRE, gp, H, JJ] = dynare_identification(options_ident, pdraws0) +% +% INPUTS +% o options_ident [structure] identification options +% o pdraws0 [matrix] optional: matrix of MC sample of model params. +% +% OUTPUTS +% o pdraws [matrix] matrix of MC sample of model params used +% o TAU, [matrix] MC sample of entries in the model solution (stacked vertically) +% o GAM, [matrix] MC sample of entries in the moments (stacked vertically) +% o LRE, [matrix] MC sample of entries in LRE model (stacked vertically) +% o gp, [matrix] derivatives of the Jacobian (LRE model) +% o H, [matrix] derivatives of the model solution +% o JJ [matrix] derivatives of the moments +% +% SPECIAL REQUIREMENTS +% None % main % @@ -27,6 +44,7 @@ else warning off MATLAB:dividebyzero end +fname_ = M_.fname; 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); @@ -44,7 +62,7 @@ if isempty(estim_params_), prior_exist=0; else prior_exist=1; -end +end iload = options_ident.load_ident_files; advanced = options_ident.advanced; @@ -59,8 +77,9 @@ options_.Schur_vec_tol = 1.e-8; options_ = set_default_option(options_,'datafile',[]); options_.mode_compute = 0; -[data,rawdata]=dynare_estimation_init([],1); -% computes a first linear solution to set up various variables +options_.plot_priors = 0; +[data,rawdata]=dynare_estimation_init([],fname_,1); + @@ -68,10 +87,12 @@ SampleSize = options_ident.prior_mc; % results = prior_sampler(0,M_,bayestopt_,options_,oo_); -if options_ident.prior_range - prior_draw(1,1); -else - prior_draw(1); +if prior_exist + if options_ident.prior_range + prior_draw(1,1); + else + prior_draw(1); + end end if ~(exist('sylvester3mr','file')==2), @@ -112,7 +133,7 @@ disp(' ') disp(['==== Identification analysis ====' ]), disp(' ') -if iload <=0, +if iload <=0, iteration = 0; burnin_iteration = 0; @@ -133,14 +154,18 @@ if iload <=0, while iteration < SampleSize, loop_indx = loop_indx+1; if prior_exist, - if nargin==2, - if burnin_iteration>=BurninSampleSize, - params = pdraws0(iteration+1,:); - else - params = pdraws0(burnin_iteration+1,:); - end + if SampleSize==1, + params = set_prior(estim_params_,M_,options_)'; else - params = prior_draw(); + if nargin==2, + if burnin_iteration>=BurninSampleSize, + params = pdraws0(iteration+1,:); + else + params = pdraws0(burnin_iteration+1,:); + end + else + params = prior_draw(); + end end set_all_parameters(params); else @@ -151,11 +176,6 @@ if iload <=0, if info(1)==0, oo0=oo_; - % [Aa,Bb] = kalman_transition_matrix(oo0.dr, ... - % bayestopt_.restrict_var_list, ... - % bayestopt_.restrict_columns, ... - % bayestopt_.restrict_aux, M_.exo_nbr); - % tau=[vec(Aa); dyn_vech(Bb*M_.Sigma_e*Bb')]; 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,1); @@ -164,7 +184,7 @@ if iload <=0, burnin_iteration = burnin_iteration + 1; pdraws(burnin_iteration,:) = params; TAU(:,burnin_iteration)=tau; - LRE(:,burnin_iteration)=[oo_.dr.ys(oo_.dr.order_var); vec(g1)]; + LRE(:,burnin_iteration)=[oo_.dr.ys(oo_.dr.order_var); vec(g1)]; [gam,stationary_vars] = th_autocovariances(oo0.dr,bayestopt_.mfys,M_,options_); if exist('OCTAVE_VERSION') warning('off', 'Octave:divide-by-zero') @@ -206,7 +226,7 @@ if iload <=0, end if iteration, - [JJ, H, gam, gp] = getJJ(A, B, M_,oo0,options_,0,indx,indexo,bayestopt_.mf2,nlags,useautocorr); + [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)); @@ -239,24 +259,24 @@ if iload <=0, % 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; - siLREmean = abs(siLRE)./SampleSize; - derJmean = (siJ)./SampleSize; - derHmean = (siH)./SampleSize; - derLREmean = (siLRE)./SampleSize; + siJmean = abs(siJ)./SampleSize; + siHmean = abs(siH)./SampleSize; + siLREmean = abs(siLRE)./SampleSize; + derJmean = (siJ)./SampleSize; + derHmean = (siH)./SampleSize; + derLREmean = (siLRE)./SampleSize; else - siJmean = abs(siJ)./SampleSize+siJmean; - siHmean = abs(siH)./SampleSize+siHmean; - siLREmean = abs(siLRE)./SampleSize+siLREmean; - derJmean = (siJ)./SampleSize+derJmean; - derHmean = (siH)./SampleSize+derHmean; - derLREmean = (siLRE)./SampleSize+derLREmean; + siJmean = abs(siJ)./SampleSize+siJmean; + siHmean = abs(siH)./SampleSize+siHmean; + siLREmean = abs(siLRE)./SampleSize+siLREmean; + derJmean = (siJ)./SampleSize+derJmean; + derHmean = (siH)./SampleSize+derHmean; + derLREmean = (siLRE)./SampleSize+derLREmean; end pdraws(iteration,:) = params; normH = max(abs(stoH(:,:,run_index))')'; @@ -265,74 +285,108 @@ if iload <=0, % normH = TAU(:,iteration); % normJ = GAM(:,iteration); % normLRE = LRE(:,iteration); - [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)), bayestopt_); - % identification_checks(H(indH,:),JJ(indJJ,:), gp(indLRE,:), bayestopt_); + [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); ide_strength_J(iteration,:)=NaN(1,nparam); - if iteration ==1, + 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]); + end cmm = simulated_moment_uncertainty(indJJ, periods, replic); + % 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.*normaliz(indok)'; + ide_strength_J(iteration,indok) = (1./(sqrt(diag(inv(MIM)))./normaliz(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)']); + normaliz(indok) = sqrt(diag(inv(MIM)))'; + % inok = find((abs(GAM(:,iteration))==0)); + % isok = find((abs(GAM(:,iteration)))); + % quant(isok,:) = siJ(isok,:)./repmat(GAM(isok,iteration),1,nparam); + % quant(inok,:) = siJ(inok,:)./repmat(mean(abs(GAM(:,iteration))),length(inok),nparam); + 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))==0)); + 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))==0)); + 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, - % 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 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', ... - 'siHmean', 'siJmean', 'siLREmean', 'derHmean', 'derJmean', 'derLREmean', 'TAU', 'GAM', 'LRE') + + + save([IdentifDirectoryName '/' M_.fname '_identif.mat'], 'pdraws', 'idemodel', 'idemoments', 'idelre', 'indJJ', 'indH', 'indLRE', ... + 'siHmean', 'siJmean', 'siLREmean', 'derHmean', 'derJmean', 'derLREmean', 'TAU', 'GAM', 'LRE') else - load([IdentifDirectoryName '/' M_.fname '_identif'], 'pdraws', 'idemodel', 'idemoments', 'idelre', 'indJJ', 'indH', 'indLRE', ... - 'siHmean', 'siJmean', 'siLREmean', 'derHmean', 'derJmean', 'derLREmean', 'TAU', 'GAM', 'LRE') + load([IdentifDirectoryName '/' M_.fname '_identif'], 'pdraws', 'idemodel', 'idemoments', 'idelre', 'indJJ', 'indH', 'indLRE', ... + 'siHmean', 'siJmean', 'siLREmean', 'derHmean', 'derJmean', 'derLREmean', 'TAU', 'GAM', 'LRE') identFiles = dir([IdentifDirectoryName '/' M_.fname '_identif_*']); options_ident.prior_mc=size(pdraws,1); - SampleSize = options_ident.prior_mc; + SampleSize = options_ident.prior_mc; options_.options_ident = options_ident; if iload>1, idemodel0=idemodel; @@ -341,7 +395,7 @@ else 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') + load([IdentifDirectoryName '/' M_.fname '_identif_' int2str(file_index)], 'stoH', 'stoJJ', 'stoLRE') for index=1:size(stoH,3), iteration = iteration+1; normH = max(abs(stoH(:,:,index))')'; @@ -350,26 +404,26 @@ else % normH = TAU(:,iteration); % normJ = GAM(:,iteration); % normLRE = LRE(:,iteration); - [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)] = ... + [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(stoH(:,:,index)./normH(:,ones(nparam,1)), ... - stoJJ(:,:,index)./normJ(:,ones(nparam,1)), ... - stoLRE(:,:,index)./normLRE(:,ones(size(stoLRE,2),1)), bayestopt_); + stoJJ(:,:,index)./normJ(:,ones(nparam,1)), ... + stoLRE(:,:,index)./normLRE(:,ones(size(stoLRE,2),1))); waitbar(iteration/SampleSize,h,['MC Identification checks ',int2str(iteration),'/',int2str(SampleSize)]) end end close(h); - save([IdentifDirectoryName '/' M_.fname '_identif'], 'idemodel', 'idemoments', 'idelre', '-append') + save([IdentifDirectoryName '/' M_.fname '_identif.mat'], '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') + 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); @@ -380,7 +434,7 @@ else 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 @@ -422,7 +476,7 @@ if SampleSize>1, if j>offset indd = 1:length(siLREmean(:,j-offset)); tstLREmean(indd,j-offset) = abs(derLREmean(indd,j-offset))./siLREmean(indd,j-offset); - end + end end end @@ -431,114 +485,18 @@ if nargout>3 && iload, filnam = dir([IdentifDirectoryName '/' M_.fname '_identif_*.mat']); H=[]; JJ = []; - gp = []; + gp = []; for j=1:length(filnam), load([IdentifDirectoryName '/' M_.fname '_identif_',int2str(j),'.mat']); H = cat(3,H, stoH(:,abs(iload),:)); JJ = cat(3,JJ, stoJJ(:,abs(iload),:)); gp = cat(3,gp, stoLRE(:,abs(iload),:)); - + end end -% mTAU = mean(TAU'); -% mGAM = mean(GAM'); -% sTAU = std(TAU'); -% sGAM = std(GAM'); -% if nargout>=3, -% GAM0=GAM; -% end -% if useautocorr, -% idiag = find(dyn_vech(eye(size(options_.varobs,1)))); -% GAM(idiag,:) = GAM(idiag,:)./(sGAM(idiag)'*ones(1,SampleSize)); -% % siJmean(idiag,:) = siJmean(idiag,:)./(sGAM(idiag)'*ones(1,nparam)); -% % siJmean = siJmean./(max(siJmean')'*ones(size(params))); -% end -% -% [pcc, dd] = eig(cov(GAM')); -% [latent, isort] = sort(-diag(dd)); -% latent = -latent; -% pcc=pcc(:,isort); -% siPCA = (siJmean'*abs(pcc')).^2'; -% siPCA = siPCA./(max(siPCA')'*ones(1,nparam)).*(latent*ones(1,nparam)); -% siPCA = sum(siPCA,1); -% siPCA = siPCA./max(siPCA); -% -% [pcc, dd] = eig(corrcoef(GAM')); -% [latent, isort] = sort(-diag(dd)); -% latent = -latent; -% pcc=pcc(:,isort); -% siPCA2 = (siJmean'*abs(pcc')).^2'; -% siPCA2 = siPCA2./(max(siPCA2')'*ones(1,nparam)).*(latent*ones(1,nparam)); -% siPCA2 = sum(siPCA2,1); -% siPCA2 = siPCA2./max(siPCA2); -% -% [pcc, dd] = eig(cov(TAU')); -% [latent, isort] = sort(-diag(dd)); -% latent = -latent; -% pcc=pcc(:,isort); -% siHPCA = (siHmean'*abs(pcc')).^2'; -% siHPCA = siHPCA./(max(siHPCA')'*ones(1,nparam)).*(latent*ones(1,nparam)); -% siHPCA = sum(siHPCA,1); -% siHPCA = siHPCA./max(siHPCA); -% -% [pcc, dd] = eig(corrcoef(TAU')); -% [latent, isort] = sort(-diag(dd)); -% latent = -latent; -% pcc=pcc(:,isort); -% siHPCA2 = (siHmean'*abs(pcc')).^2'; -% siHPCA2 = siHPCA2./(max(siHPCA2')'*ones(1,nparam)).*(latent*ones(1,nparam)); -% siHPCA2 = sum(siHPCA2,1); -% siHPCA2 = siHPCA2./max(siHPCA2); - - disp_identification(pdraws, idemodel, idemoments, name, advanced) -% figure, -% % myboxplot(siPCA(1:(max(find(cumsum(latent)./length(indJJ)<0.99))+1),:)) -% subplot(221) -% bar(siHPCA) -% % set(gca,'ylim',[0 1]) -% set(gca,'xticklabel','') -% set(gca,'xlim',[0.5 nparam+0.5]) -% for ip=1:nparam, -% text(ip,-0.02,name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none') -% end -% title('Sensitivity in TAU''s PCA') -% -% subplot(222) -% % myboxplot(siPCA(1:(max(find(cumsum(latent)./length(indJJ)<0.99))+1),:)) -% bar(siHPCA2) -% % set(gca,'ylim',[0 1]) -% set(gca,'xticklabel','') -% set(gca,'xlim',[0.5 nparam+0.5]) -% for ip=1:nparam, -% text(ip,-0.02,name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none') -% end -% title('Sensitivity in standardized TAU''s PCA') -% -% -% subplot(223) -% % myboxplot(siPCA(1:(max(find(cumsum(latent)./length(indJJ)<0.99))+1),:)) -% bar(siPCA) -% % set(gca,'ylim',[0 1]) -% set(gca,'xticklabel','') -% set(gca,'xlim',[0.5 nparam+0.5]) -% for ip=1:nparam, -% text(ip,-0.02,name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none') -% end -% title('Sensitivity in moments'' PCA') -% -% subplot(224) -% % myboxplot(siPCA(1:(max(find(cumsum(latent)./length(indJJ)<0.99))+1),:)) -% bar(siPCA2) -% % set(gca,'ylim',[0 1]) -% set(gca,'xticklabel','') -% set(gca,'xlim',[0.5 nparam+0.5]) -% for ip=1:nparam, -% 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'), @@ -563,7 +521,7 @@ if advanced, 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'); @@ -586,10 +544,10 @@ if advanced, eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_LRE']); 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 @@ -609,7 +567,7 @@ if advanced, % 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); % [ss, is] = sort(mmm); @@ -624,7 +582,7 @@ if advanced, myboxplot(log(siHnorm(:,is))) else bar(siHnorm(:,is)) - end + end % set(gca,'ylim',[0 1.05]) set(gca,'xticklabel','') dy = get(gca,'ylim'); @@ -633,7 +591,7 @@ if advanced, 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'); @@ -672,15 +630,16 @@ end if SampleSize>1, myboxplot(log(ide_strength_J(:,is))) else - bar(ide_strength_J(:,is)) -end + bar(log(ide_strength_J(:,is))) +end % set(gca,'ylim',[0 1.05]) +set(gca,'xlim',[0 nparam+1]) 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 +% for ip=1:nparam, +% text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none') +% end title('Identification strength in the moments') subplot(212) @@ -699,6 +658,7 @@ else bar(siJnorm(:,is)) end % set(gca,'ylim',[0 1.05]) +set(gca,'xlim',[0 nparam+1]) set(gca,'xticklabel','') dy = get(gca,'ylim'); % dy=dy(2)-dy(1); @@ -730,30 +690,55 @@ title('Sensitivity in the moments') % eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_moments']); % if options_.nograph, close(gcf); end -if SampleSize==1, +if SampleSize==1 && advanced, % identificaton patterns [U,S,V]=svd(siJ./normJ(:,ones(nparam,1)),0); - figure('name','Identification patterns (moments)'), + 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'); + end for j=1:min(nparam,8), - subplot(2,4,j), + 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); - if j==1, ylabel('SMALLEST singular values'), end, +% if j==4 || j==nparam, +% xlabel('SMALLEST singular values'), +% end, else bar(abs(V(:,j))), Stit = S(j,j); - if j==5, ylabel('LARGEST singular values'), end, +% if j==8 || j==nparam, +% xlabel('LARGEST singular values'), +% end, 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 - saveas(gcf,[IdentifDirectoryName,'/',M_.fname,'_ident_pattern']) - eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_pattern']); - eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_pattern']); + 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 if advanced, @@ -773,17 +758,17 @@ if advanced, eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_COND']); if options_.nograph, close(gcf); end end - - + + pco = NaN(np,np); - for j=1:np, + for j=1:np, if SampleSize>1 - pco(j+1:end,j) = mean(abs(squeeze(idelre.Pco(j+1:end,j,:))')); + pco(j+1:end,j) = mean(abs(squeeze(idelre.Pco(j+1:end,j,:))')); else - pco(j+1:end,j) = abs(idelre.Pco(j+1:end,j))'; + pco(j+1:end,j) = abs(idelre.Pco(j+1:end,j))'; end end - figure('name','Pairwise correlations in the LRE model'), + figure('name','Pairwise correlations in the LRE model'), imagesc(pco',[0 1]); set(gca,'xticklabel','') set(gca,'yticklabel','') @@ -796,16 +781,16 @@ if advanced, 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, + for j=1:nparam, if SampleSize>1 - pco(j+1:end,j) = mean(abs(squeeze(idemodel.Pco(j+1:end,j,:))')); + 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))'; + pco(j+1:end,j) = abs(idemodel.Pco(j+1:end,j))'; end end - figure('name','Pairwise correlations in the model'), + figure('name','Pairwise correlations in the model'), imagesc(pco',[0 1]); set(gca,'xticklabel','') set(gca,'yticklabel','') @@ -818,15 +803,15 @@ if advanced, 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, + + for j=1:nparam, if SampleSize>1 - pco(j+1:end,j) = mean(abs(squeeze(idemoments.Pco(j+1:end,j,:))')); + 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))'; + pco(j+1:end,j) = abs(idemoments.Pco(j+1:end,j))'; end end - figure('name','Pairwise correlations in the moments'), + figure('name','Pairwise correlations in the moments'), imagesc(pco',[0 1]); set(gca,'xticklabel','') set(gca,'yticklabel','') @@ -840,110 +825,6 @@ if advanced, 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') diff --git a/matlab/identification_checks.m b/matlab/identification_checks.m index a8b60b3bb..02c55a2cc 100644 --- a/matlab/identification_checks.m +++ b/matlab/identification_checks.m @@ -1,4 +1,38 @@ -function [McoH, McoJ, McoGP, PcoH, PcoJ, PcoGP, condH, condJ, condGP, eH, eJ, eGP, ind01, ind02, indnoH, indnoJ, ixnoH, ixnoJ] = identification_checks(H,JJ, gp, bayestopt_) +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) +% 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 +% +% 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 zero columns of H +% o ind02 [array] binary indicator for 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 +% +% SPECIAL REQUIREMENTS +% None % Copyright (C) 2008-2011 Dynare Team % @@ -20,54 +54,44 @@ 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 at theta +% 1. check rank of H, JJ, gp at theta npar = size(H,2); -npar0 = size(gp,2); +npar0 = size(gp,2); % shocks do not enter jacobian indnoH = zeros(1,npar); indnoJ = zeros(1,npar); indnoLRE = zeros(1,npar0); -ind1 = find(vnorm(H)>=eps); + +% H matrix +ind1 = find(vnorm(H)>=eps); % take non-zero columns H1 = H(:,ind1); -covH = H1'*H1; -sdH = sqrt(diag(covH)); -sdH = sdH*sdH'; -% [e1,e2] = eig( (H1'*H1)./sdH ); [eu,e2,e1] = svd( H1, 0 ); eH = zeros(npar,npar); % eH(ind1,:) = e1; -eH(ind1,length(find(vnorm(H)==0))+1:end) = e1; -eH(find(vnorm(H)==0),1:length(find(vnorm(H)==0)))=eye(length(find(vnorm(H)==0))); +eH(ind1,length(find(vnorm(H)=eps); +ind2 = find(vnorm(JJ)>=eps); % take non-zero columns JJ1 = JJ(:,ind2); -covJJ = JJ1'*JJ1; -sdJJ = sqrt(diag(covJJ)); -sdJJ = sdJJ*sdJJ'; -% [ee1,ee2] = eig( (JJ1'*JJ1)./sdJJ ); [eu,ee2,ee1] = svd( JJ1, 0 ); -% eJ = NaN(npar,length(ind2)); eJ = zeros(npar,npar); -eJ(ind2,length(find(vnorm(JJ)==0))+1:end) = ee1; -eJ(find(vnorm(JJ)==0),1:length(find(vnorm(JJ)==0)))=eye(length(find(vnorm(JJ)==0))); -condJJ = cond(JJ1'*JJ1); +eJ(ind2,length(find(vnorm(JJ)=eps); +ind3 = find(vnorm(gp)>=eps); % take non-zero columns gp1 = gp(:,ind3); covgp = gp1'*gp1; sdgp = sqrt(diag(covgp)); sdgp = sdgp*sdgp'; -[ex1,ex2] = eig( (gp1'*gp1)./sdgp ); -% eJ = NaN(npar,length(ind2)); +[eu,ex2,ex1] = svd(gp1, 0 ); eGP = zeros(npar0,npar0); -eGP(ind3,length(find(vnorm(gp)==0))+1:end) = ex1; -eGP(find(vnorm(gp)==0),1:length(find(vnorm(gp)==0)))=eye(length(find(vnorm(gp)==0))); +eGP(ind3,length(find(vnorm(gp) 1.e-6 )}; - indnoH(ixno,:) = (abs(e1(:,e0(j))) > 1.e-6 )'; - % disp('Perfectly collinear parameters') - % disp(bayestopt_.name(indnoH{ixno})) - % disp(' ') - % ind01(indnoH{ixno})=0; + indnoH(ixno,ind1) = (abs(e1(:,e0(j))) > 1.e-6 )'; end else % rank(H)==length(theta), go to 2 % 2. check rank of J @@ -135,22 +146,15 @@ if rankJ 1.e-6) )}; - indnoJ(ixno,:) = (abs(ee1(:,ee0(j))) > 1.e-6)'; - % disp('Perfectly collinear parameters in moments J') - % disp(bayestopt_.name(indnoJ{ixno})) - % disp(' ') - % ind02(indnoJ{ixno})=0; + indnoJ(ixno,ind2) = (abs(ee1(:,ee0(j))) > 1.e-6)'; end else %rank(J)==length(theta) => % disp('All parameters are identified at theta by the moments included in J') @@ -189,11 +193,6 @@ for ii = 1:size(gp1,2); end -% ind01 = zeros(npar,1); -% ind02 = zeros(npar,1); -% ind01(ind1) = 1; -% ind02(ind2) = 1; -