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

🐛 identification do plots and display results only for nonempty objects 

fixes #1860
mr#2067
Willi Mutschler 2022-07-26 10:22:45 +02:00
parent 8c01912a50
commit 30a6d35f5a
No known key found for this signature in database
GPG Key ID: 91E724BF17A73F6D
7 changed files with 122 additions and 64 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
matlab/disp_identification.m
View File

@ -17,6 +17,10 @@ function disp_identification(pdraws, ide_reducedform, ide_moments, ide_spectrum,
% (Komunjer and Ng, 2011).
% name: [totparam_nbr by 1] string cell of parameter names
% options_ident: [structure] identification options
% error_indicator [structure]
% boolean information on errors (1 is an error, 0 is no error)
% while computing the criteria
% -------------------------------------------------------------------------
% OUTPUTS
% * all output is printed on the command line
@ -97,7 +101,7 @@ for jide = 1:4
elseif options_ident.order == 3
strMeaning = 'Jacobian of first-order minimal system and third-order accurate mean';
end
if ~no_identification_minimal && ~(length(ide_minimal.minimal_state_space)==1 && ide_minimal.minimal_state_space==0)
if ~no_identification_minimal && ~(length(ide_minimal.minimal_state_space)==1 && ide_minimal.minimal_state_space==0) && isfield(ide_minimal,'dMINIMAL')
noidentification = 0; ide = ide_minimal;
if SampleSize == 1
Jacob = ide.dMINIMAL;
@ -261,7 +265,7 @@ end
%% Advanced identification patterns
if SampleSize==1 && options_ident.advanced
if SampleSize==1 && options_ident.advanced && ~isempty(fieldnames(ide_moments))
skipline()
for j=1:size(ide_moments.cosndMOMENTS,2)
pax=NaN(totparam_nbr,totparam_nbr);

36
matlab/dynare_identification.m
View File

@ -481,7 +481,7 @@ if iload <=0
end
options_ident.tittxt = parameters; %title text for graphs and figures
% perform identification analysis for single point
[ide_moments_point, ide_spectrum_point, ide_minimal_point, ide_hess_point, ide_reducedform_point, ide_dynamic_point, derivatives_info_point, info] = ...
[ide_moments_point, ide_spectrum_point, ide_minimal_point, ide_hess_point, ide_reducedform_point, ide_dynamic_point, derivatives_info_point, info, error_indicator_point] = ...
identification_analysis(params, indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1); %the 1 at the end implies initialization of persistent variables
if info(1)~=0
% there are errors in the solution algorithm
@ -498,7 +498,7 @@ if iload <=0
params = prior_draw();
options_ident.tittxt = 'Random_prior_params'; %title text for graphs and figures
% perform identification analysis
[ide_moments_point, ide_spectrum_point, ide_minimal_point, ide_hess_point, ide_reducedform_point, ide_dynamic_point, derivatives_info_point, info] = ...
[ide_moments_point, ide_spectrum_point, ide_minimal_point, ide_hess_point, ide_reducedform_point, ide_dynamic_point, derivatives_info_point, info, error_indicator_point] = ...
identification_analysis(params, indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1);
end
end
@ -525,7 +525,7 @@ if iload <=0
save([IdentifDirectoryName '/' fname '_' parameters '_identif.mat'], 'ide_moments_point', 'ide_spectrum_point', 'ide_minimal_point', 'ide_hess_point', 'ide_reducedform_point', 'ide_dynamic_point', 'store_options_ident');
% display results of identification analysis
disp_identification(params, ide_reducedform_point, ide_moments_point, ide_spectrum_point, ide_minimal_point, name, options_ident);
if ~options_ident.no_identification_strength && ~options_.nograph
if ~options_ident.no_identification_strength && ~options_.nograph && ~error_indicator_point.identification_strength && ~error_indicator_point.identification_moments
% plot (i) identification strength and sensitivity measure based on the moment information matrix and (ii) plot advanced analysis graphs
plot_identification(params, ide_moments_point, ide_hess_point, ide_reducedform_point, ide_dynamic_point, options_ident.advanced, parameters, name, IdentifDirectoryName, parameters_TeX, name_tex);
end
@ -589,7 +589,7 @@ if iload <=0
end
% preallocate storage for moments
if ~options_MC.no_identification_moments
if ~options_MC.no_identification_moments && ~isempty(fieldnames(ide_moments))
STO_si_dMOMENTS = zeros([size(ide_moments.si_dMOMENTS, 1), totparam_nbr, MAX_RUNS_BEFORE_SAVE_TO_FILE]);
STO_MOMENTS = zeros(size(ide_moments.MOMENTS, 1), SampleSize);
IDE_MOMENTS.ind_dMOMENTS = ide_moments.ind_dMOMENTS;
@ -608,7 +608,7 @@ if iload <=0
end
% preallocate storage for spectrum
if ~options_MC.no_identification_spectrum
if ~options_MC.no_identification_spectrum && ~isempty(fieldnames(ide_spectrum))
STO_dSPECTRUM = zeros([size(ide_spectrum.dSPECTRUM, 1), size(ide_spectrum.dSPECTRUM, 2), MAX_RUNS_BEFORE_SAVE_TO_FILE]);
IDE_SPECTRUM.ind_dSPECTRUM = ide_spectrum.ind_dSPECTRUM;
IDE_SPECTRUM.ino = zeros(SampleSize, 1);
@ -623,7 +623,7 @@ if iload <=0
end
% preallocate storage for minimal system
if ~options_MC.no_identification_minimal
if ~options_MC.no_identification_minimal && ~isempty(fieldnames(ide_minimal)) && ide_minimal.minimal_state_space==1
STO_dMINIMAL = zeros([size(ide_minimal.dMINIMAL, 1), size(ide_minimal.dMINIMAL, 2), MAX_RUNS_BEFORE_SAVE_TO_FILE]);
IDE_MINIMAL.ind_dMINIMAL = ide_minimal.ind_dMINIMAL;
IDE_MINIMAL.ino = zeros(SampleSize, 1);
@ -801,21 +801,21 @@ if iload <=0
% note that this is not the same si_dDYNAMICnorm as computed in identification_analysis
% given that we have the MC sample of the Jacobians, we also normalize by the std of the sample of Jacobian entries, to get a fully standardized sensitivity measure
si_dDYNAMICnorm(iter,:) = vnorm(STO_si_dDYNAMIC(:,:,irun)./repmat(normalize_STO_DYNAMIC,1,totparam_nbr-(stderrparam_nbr+corrparam_nbr))).*normaliz1((stderrparam_nbr+corrparam_nbr)+1:end);
if ~options_MC.no_identification_reducedform
if ~options_MC.no_identification_reducedform && ~isempty(STO_si_dREDUCEDFORM)
% note that this is not the same si_dREDUCEDFORMnorm as computed in identification_analysis
% given that we have the MC sample of the Jacobians, we also normalize by the std of the sample of Jacobian entries, to get a fully standardized sensitivity measure
si_dREDUCEDFORMnorm(iter,:) = vnorm(STO_si_dREDUCEDFORM(:,:,irun)./repmat(normalize_STO_REDUCEDFORM,1,totparam_nbr)).*normaliz1;
end
if ~options_MC.no_identification_moments
if ~options_MC.no_identification_moments && ~isempty(STO_si_dMOMENTS)
% note that this is not the same si_dMOMENTSnorm as computed in identification_analysis
% given that we have the MC sample of the Jacobians, we also normalize by the std of the sample of Jacobian entries, to get a fully standardized sensitivity measure
si_dMOMENTSnorm(iter,:) = vnorm(STO_si_dMOMENTS(:,:,irun)./repmat(normalize_STO_MOMENTS,1,totparam_nbr)).*normaliz1;
end
if ~options_MC.no_identification_spectrum
if ~options_MC.no_identification_spectrum && ~isempty(STO_dSPECTRUM)
% note that this is not the same dSPECTRUMnorm as computed in identification_analysis
dSPECTRUMnorm(iter,:) = vnorm(STO_dSPECTRUM(:,:,irun)); %not yet used
end
if ~options_MC.no_identification_minimal
if ~options_MC.no_identification_minimal && ~isempty(STO_dMINIMAL)
% note that this is not the same dMINIMALnorm as computed in identification_analysis
dMINIMALnorm(iter,:) = vnorm(STO_dMINIMAL(:,:,irun)); %not yet used
end
@ -874,7 +874,7 @@ if iload
%if previous analysis is loaded
fprintf(['Testing %s\n',parameters]);
disp_identification(ide_hess_point.params, ide_reducedform_point, ide_moments_point, ide_spectrum_point, ide_minimal_point, name, options_ident);
if ~options_.nograph
if ~options_.nograph && ~error_indicator_point.identification_strength && ~error_indicator_point.identification_moments
% plot (i) identification strength and sensitivity measure based on the sample information matrix and (ii) advanced analysis graphs
plot_identification(ide_hess_point.params, ide_moments_point, ide_hess_point, ide_reducedform_point, ide_dynamic_point, options_ident.advanced, parameters, name, IdentifDirectoryName, [], name_tex);
end
@ -888,7 +888,7 @@ if SampleSize > 1
options_ident.advanced = 0;
disp_identification(pdraws, IDE_REDUCEDFORM, IDE_MOMENTS, IDE_SPECTRUM, IDE_MINIMAL, name, options_ident);
options_ident.advanced = advanced0; % reset advanced setting
if ~options_.nograph
if ~options_.nograph && isfield(ide_hess_point,'ide_strength_dMOMENTS')
% plot (i) identification strength and sensitivity measure based on the sample information matrix and (ii) advanced analysis graphs
plot_identification(pdraws, IDE_MOMENTS, ide_hess_point, IDE_REDUCEDFORM, IDE_DYNAMIC, options_ident.advanced, 'MC sample ', name, IdentifDirectoryName, [], name_tex);
end
@ -906,14 +906,14 @@ if SampleSize > 1
fprintf('\nTesting %s.\n',tittxt);
if ~iload
options_ident.tittxt = tittxt; %title text for graphs and figures
[ide_moments_max, ide_spectrum_max, ide_minimal_max, ide_hess_max, ide_reducedform_max, ide_dynamic_max, derivatives_info_max, info_max] = ...
[ide_moments_max, ide_spectrum_max, ide_minimal_max, ide_hess_max, ide_reducedform_max, ide_dynamic_max, derivatives_info_max, info_max, error_indicator_max] = ...
identification_analysis(pdraws(jmax,:), indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1); %the 1 at the end initializes some persistent variables
save([IdentifDirectoryName '/' fname '_identif.mat'], 'ide_hess_max', 'ide_moments_max', 'ide_spectrum_max', 'ide_minimal_max','ide_reducedform_max', 'ide_dynamic_max', 'jmax', '-append');
end
advanced0 = options_ident.advanced; options_ident.advanced = 1; % make sure advanced setting is on
disp_identification(pdraws(jmax,:), ide_reducedform_max, ide_moments_max, ide_spectrum_max, ide_minimal_max, name, options_ident);
options_ident.advanced = advanced0; %reset advanced setting
if ~options_.nograph
if ~options_.nograph && ~error_indicator_max.identification_strength && ~error_indicator_max.identification_moments
% plot (i) identification strength and sensitivity measure based on the sample information matrix and (ii) advanced analysis graphs
plot_identification(pdraws(jmax,:), ide_moments_max, ide_hess_max, ide_reducedform_max, ide_dynamic_max, 1, tittxt, name, IdentifDirectoryName, tittxt, name_tex);
end
@ -924,14 +924,14 @@ if SampleSize > 1
fprintf('Testing %s.\n',tittxt);
if ~iload
options_ident.tittxt = tittxt; %title text for graphs and figures
[ide_moments_min, ide_spectrum_min, ide_minimal_min, ide_hess_min, ide_reducedform_min, ide_dynamic_min, derivatives_info_min, info_min] = ...
[ide_moments_min, ide_spectrum_min, ide_minimal_min, ide_hess_min, ide_reducedform_min, ide_dynamic_min, derivatives_info_min, info_min, error_indicator_min] = ...
identification_analysis(pdraws(jmin,:), indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1); %the 1 at the end initializes persistent variables
save([IdentifDirectoryName '/' fname '_identif.mat'], 'ide_hess_min', 'ide_moments_min','ide_spectrum_min','ide_minimal_min','ide_reducedform_min', 'ide_dynamic_min', 'jmin', '-append');
end
advanced0 = options_ident.advanced; options_ident.advanced = 1; % make sure advanced setting is on
disp_identification(pdraws(jmin,:), ide_reducedform_min, ide_moments_min, ide_spectrum_min, ide_minimal_min, name, options_ident);
options_ident.advanced = advanced0; %reset advanced setting
if ~options_.nograph
if ~options_.nograph && ~error_indicator_min.identification_strength && ~error_indicator_min.identification_moments
% plot (i) identification strength and sensitivity measure based on the sample information matrix and (ii) advanced analysis graphs
plot_identification(pdraws(jmin,:),ide_moments_min,ide_hess_min,ide_reducedform_min,ide_dynamic_min,1,tittxt,name,IdentifDirectoryName,tittxt,name_tex);
end
@ -946,13 +946,13 @@ if SampleSize > 1
fprintf('\nTesting %s.\n',tittxt);
if ~iload
options_ident.tittxt = tittxt; %title text for graphs and figures
[ide_moments_(j), ide_spectrum_(j), ide_minimal_(j), ide_hess_(j), ide_reducedform_(j), ide_dynamic_(j), derivatives_info_(j), info_resolve] = ...
[ide_moments_(j), ide_spectrum_(j), ide_minimal_(j), ide_hess_(j), ide_reducedform_(j), ide_dynamic_(j), derivatives_info_(j), info_resolve, error_indicator_j] = ...
identification_analysis(pdraws(jcrit(j),:), indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1);
end
advanced0 = options_ident.advanced; options_ident.advanced = 1; %make sure advanced setting is on
disp_identification(pdraws(jcrit(j),:), ide_reducedform_(j), ide_moments_(j), ide_spectrum_(j), ide_minimal_(j), name, options_ident);
options_ident.advanced = advanced0; % reset advanced
if ~options_.nograph
if ~options_.nograph && ~error_indicator_j.identification_strength && ~error_indicator_j.identification_moments
% plot (i) identification strength and sensitivity measure based on the sample information matrix and (ii) advanced analysis graphs
plot_identification(pdraws(jcrit(j),:), ide_moments_(j), ide_hess_(j), ide_reducedform_(j), ide_dynamic_(j), 1, tittxt, name, IdentifDirectoryName, tittxt, name_tex);
end

10
matlab/identification_analysis.m
View File

@ -29,6 +29,8 @@ function [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide
% 0: prior does not exist. Identification is checked for all stderr and model parameters, but no corr parameters
% * init [integer]
% flag for initialization of persistent vars. This is needed if one may want to make more calls to identification in the same mod file
% * error_indicator [structure]
% boolean information on errors (1 is an error, 0 is no error) while computing the criteria stored in fields identification_strength, identification_reducedform, identification_moments, identification_minimal, identification_spectrum
% -------------------------------------------------------------------------
% OUTPUTS
% * ide_moments [structure]
@ -209,7 +211,7 @@ if info(1) == 0 %no errors in solution
warning_MOMENTS = [warning_MOMENTS ' The number of moments with non-zero derivative is smaller than the number of parameters up to 10 lags.\n'];
warning_MOMENTS = [warning_MOMENTS ' Either reduce the list of parameters, or further increase ar, or increase number of observables.\n'];
warning_MOMENTS = [warning_MOMENTS ' Skip identification analysis based on moments.\n'];
warning_MOMENTS = [warning_MOMENTS ' Skip identification strenght analysis.\n'];
warning_MOMENTS = [warning_MOMENTS ' Skip identification strength analysis.\n'];
fprintf(warning_MOMENTS);
%set indicator to neither display nor plot dMOMENTS anymore
error_indicator.identification_moments = 1;
@ -472,7 +474,11 @@ if info(1) == 0 %no errors in solution
%here we focus on the unnormalized S and V, which is then used in plot_identification.m and for prior_mc > 1
[~, S, V] = svd(dMOMENTS(ind_dMOMENTS,:),0);
S = diag(S);
if size(S,1) == 1
S = S(1); % edge case that S is not a matrix but a row vector
else
S = diag(S);
end
S = [S;zeros(size(dMOMENTS,2)-length(ind_dMOMENTS),1)];
if totparam_nbr > 8
ide_moments.S = S([1:4, end-3:end]);

4
tests/Makefile.am
View File

@ -269,7 +269,8 @@ MODFILES = \
identification/rbc_ident/rbc_ident_std_as_structural_par.mod \
identification/rbc_ident/rbc_ident_varexo_only.mod \
identification/correlated_errors/fs2000_corr.mod \
identification/forward_looking/forward_looking.mod \
identification/forward_looking/forward_looking_empty_ghx.mod \
identification/forward_looking/forward_looking_varobs_x.mod \
simul/example1.mod \
simul/Solow_no_varexo.mod \
simul/simul_ZLB_purely_forward.mod \
@ -1330,6 +1331,7 @@ EXTRA_DIST = \
kalman/likelihood_from_dynare/fs2000ns_model.inc \
kalman/likelihood_from_dynare/fs2000ns_estimation_check.inc \
identification/as2007/G_QT.mat \
identification/forward_looking/forward_looking_common.inc \
estimation/fsdat_simul.m \
estimation/heteroskedastic_shocks/fs2000_het_model.inc \
estimation/heteroskedastic_shocks/fs2000_het_check.inc \

68
tests/identification/forward_looking/forward_looking.mod → tests/identification/forward_looking/forward_looking_common.inc Executable file → Normal file
View File

@ -1,7 +1,7 @@
% Forward-looking example model from Koop, Pesaran, Smith (2013, JBES)
% created by Willi Mutschler (@wmutschl, willi@mutschler.eu)
% =========================================================================
% Copyright © 2020 Dynare Team
% Copyright © 2020-2022 Dynare Team
%
% This file is part of Dynare.
%
@ -18,53 +18,39 @@
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% =========================================================================
var r x p;
varexo e_M e_D e_S;
varobs r x p;
var
r ${r}$ (long_name='interest rate')
x ${x}$ (long_name='output-gap')
p ${\pi}$ (long_name='inflation')
;
varexo
e_M ${\varepsilon^M}$ (long_name='monetary policy shock')
e_D ${\varepsilon^D}$ (long_name='demand shock')
e_S ${\varepsilon^S}$ (long_name='supply shock')
;
parameters
PSI ${\psi}$ (long_name='inflation elasticity Taylor rule')
TAU ${\tau}$ (long_name='intertemporal elasticity of subsitution')
BETA ${\beta}$ (long_name='discount factor')
KAPPA ${\kappa}$ (long_name='slope Phillips curve')
;
parameters PSI TAU BETA KAPPA;
PSI=1.1;
TAU=2;
BETA=0.9;
KAPPA=0.6;
PSI = 1.1;
TAU = 2;
BETA = 0.9;
KAPPA = 0.6;
model;
[name='Taylor rule']
r = PSI*p + e_M;
[name='New Keynesian IS curve']
x = x(+1) - 1/TAU*(r-p(+1)) + e_D;
[name='New Keynesian Phillips curve']
p = BETA*p(+1) + KAPPA*x + e_S;
end;
shocks;
var e_M = 1;
var e_D = 1;
var e_S = 1;
var e_M = 1;
var e_D = 1;
var e_S = 1;
end;
steady;
check;
estimated_params;
PSI, 1.1;
TAU, 2;
BETA, 0.9;
KAPPA, 0.6;
end;
identification; %this triggers sylvester3a with empty ghx
% as a side note, we have the true solution:
% [r;x;p] = TRUE_SOLUTION*[e_M;e_D;e_S] (ghx is empty)
A = [1 0 -PSI; 1/TAU 1 0; 0 -KAPPA 1];
TRUE_SOLUTION1 = inv(A);
TRUE_SOLUTION2 = 1/(KAPPA*PSI/TAU +1)*[1 KAPPA*PSI PSI;
-1/TAU 1 -PSI/TAU;
-KAPPA/TAU KAPPA 1];
% note that BETA drops out from the solution
if max(max(abs(TRUE_SOLUTION1 - oo_.dr.ghu))) > 1e-15
error('Something wrong with perturbation');
end
if max(max(abs(TRUE_SOLUTION2 - oo_.dr.ghu))) > 1e-15
error('Something wrong with perturbation');
end

35
tests/identification/forward_looking/forward_looking_empty_ghx.mod Executable file
View File

@ -0,0 +1,35 @@
@#include "forward_looking_common.inc"
PSI=1.1;
TAU=2;
BETA=0.9;
KAPPA=0.6;
steady;
check;
estimated_params;
PSI, 1.1;
TAU, 2;
BETA, 0.9;
KAPPA, 0.6;
end;
varobs r x p;
identification; %this triggers sylvester3a with empty ghx
% as a side note, we have the true solution:
% [r;x;p] = TRUE_SOLUTION*[e_M;e_D;e_S] (ghx is empty)
A = [1 0 -PSI; 1/TAU 1 0; 0 -KAPPA 1];
TRUE_SOLUTION1 = inv(A);
TRUE_SOLUTION2 = 1/(KAPPA*PSI/TAU +1)*[1 KAPPA*PSI PSI;
-1/TAU 1 -PSI/TAU;
-KAPPA/TAU KAPPA 1];
% note that BETA drops out from the solution
if max(max(abs(TRUE_SOLUTION1 - oo_.dr.ghu))) > 1e-15
error('Something wrong with perturbation');
end
if max(max(abs(TRUE_SOLUTION2 - oo_.dr.ghu))) > 1e-15
error('Something wrong with perturbation');
end

25
tests/identification/forward_looking/forward_looking_varobs_x.mod Normal file
View File

@ -0,0 +1,25 @@
@#include "forward_looking_common.inc"
PSI=1.1;
TAU=2;
BETA=0.9;
KAPPA=0.6;
steady;
check;
estimated_params;
PSI, 1.1, NORMAL_PDF, 1.1, 0.01;
TAU, 2, NORMAL_PDF, 2, 0.01;
BETA, 0.9, NORMAL_PDF, 0.9, 0.01;
KAPPA, 0.6, NORMAL_PDF, 0.6, 0.01;
end;
% with this calibration and just a single observable the identification tests run into several edge case
% due to the fact that the order condition for dMoments fails and
% the minimal system cannot be computed
varobs x;
identification;
identification(advanced=1);
identification(prior_mc=10);
identification(prior_mc=10,advanced=1);