Skip to content
GitLab
Verified Commit c4f7c416 authored by Willi Mutschler's avatar Willi Mutschler
Browse files

🐛 Fix #1694 by robust rank tolerance and histc

parent d3e90a8d
Hide whitespace changes
Inline Side-by-side
matlab/disp_identification.m
View file @ c4f7c416
......@@ -71,7 +71,7 @@ if options_ident.normalize_jacobians == 1
else
fprintf(' Normalize Jacobians: No\n');
end
fprintf(' Tolerance level for rank computations: %.0d\n',options_ident.tol_rank);
fprintf(' Tolerance level for rank computations: %s\n',num2str(options_ident.tol_rank));
fprintf(' Tolerance level for selecting nonzero columns: %.0d\n',options_ident.tol_deriv);
fprintf(' Tolerance level for selecting nonzero singular values: %.0d\n',options_ident.tol_sv);
......
matlab/dynare_identification.m
View file @ c4f7c416
......@@ -89,6 +89,7 @@ else
warning('off','MATLAB:specgraph:private:specgraph:UsingOnlyRealComponentOfComplexData');
warning('off','MATLAB:handle_graphics:exceptions:SceneNode');
warning('off','MATLAB:divideByZero');
warning('off','MATLAB:log:logOfZero');
end
%% Set all options and create objects
......@@ -136,7 +137,7 @@ options_ident = set_default_option(options_ident,'grid_nbr',5000);
error('IDENTIFICATION: You need to set an even value for ''grid_nbr''');
end
end
options_ident = set_default_option(options_ident,'tol_rank',1.e-10);
options_ident = set_default_option(options_ident,'tol_rank','robust');
% tolerance level used for rank computations
options_ident = set_default_option(options_ident,'tol_deriv',1.e-8);
% tolerance level for selecting columns of non-zero derivatives
......
matlab/identification_analysis.m
View file @ c4f7c416
......@@ -231,14 +231,14 @@ if info(1) == 0 %no errors in solution
end
ind_dDYNAMIC = (find(max(abs(dDYNAMIC'),[],1) > tol_deriv)); %index with non-zero rows
end
DYNAMIC = DYNAMIC(ind_dDYNAMIC); %focus only on non-zero entries
si_dDYNAMIC = (dDYNAMIC(ind_dDYNAMIC,:)); %focus only on non-zero rows
if ~no_identification_reducedform
REDUCEDFORM = REDUCEDFORM(ind_dREDUCEDFORM); %focus only on non-zero entries
si_dREDUCEDFORM = (dREDUCEDFORM(ind_dREDUCEDFORM,:)); %focus only on non-zero rows
end
if ~no_identification_moments
MOMENTS = MOMENTS(ind_dMOMENTS); %focus only on non-zero entries
si_dMOMENTS = (dMOMENTS(ind_dMOMENTS,:)); %focus only on non-zero derivatives
......@@ -268,6 +268,9 @@ if info(1) == 0 %no errors in solution
try
%try to compute asymptotic Hessian for identification strength analysis based on moments
if options_.order > 1
error('IDENTIFICATION STRENGTH: Analytic computation of Hessian is not available for ''order>1''. Identification strength is based on simulated moment uncertainty');
end
% reset some options for faster computations
options_.irf = 0;
options_.noprint = 1;
......@@ -282,10 +285,6 @@ if info(1) == 0 %no errors in solution
dataset_ = dseries(oo_.endo_simul(options_.varobs_id,100+1:end)',dates('1Q1'), options_.varobs); %get information on moments
derivatives_info.no_DLIK = 1;
bounds = prior_bounds(bayestopt_, options_.prior_trunc); %reset bounds as lb and ub must only be operational during mode-finding
if options_.order > 1
fprintf('IDENTIFICATION STRENGTH: Analytic computation of Hessian is not available for ''order>1''\n')
fprintf(' Identification strength is based on simulated moment uncertainty\n');
end
%note that for order>1 we do not provide any information on DT,DYss,DOM in derivatives_info, such that dsge_likelihood creates an error. Therefore the computation will be based on simulated_moment_uncertainty for order>1.
[fval, info, cost_flag, DLIK, AHess, ys, trend_coeff, M_, options_, bayestopt_, oo_] = dsge_likelihood(params', dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds, oo_, derivatives_info); %non-used output variables need to be set for octave for some reason
%note that for the order of parameters in AHess we have: stderr parameters come first, corr parameters second, model parameters third. the order within these blocks corresponds to the order specified in the estimated_params block
......@@ -315,6 +314,13 @@ if info(1) == 0 %no errors in solution
flag_score = 1; %this is used for the title in plot_identification.m
catch
%Asymptotic Hessian via simulation
if options_.order > 1
% reset some options for faster computations
options_.irf = 0;
options_.noprint = 1;
options_.SpectralDensity.trigger = 0;
options_.periods = periods+100;
end
replic = max([replic, length(ind_dMOMENTS)*3]);
cmm = simulated_moment_uncertainty(ind_dMOMENTS, periods, replic,options_,M_,oo_); %covariance matrix of moments
sd = sqrt(diag(cmm));
......
matlab/identification_checks.m
View file @ c4f7c416
......@@ -51,8 +51,8 @@ function [condX, rankX, ind0, indno, ixno, Mco, Pco, jweak, jweak_pair] = identi
if issparse(X)
X = full(X);
end
if nargin < 3 || isempty(tol_rank)
tol_rank = 1.e-10; %tolerance level used for rank computations
if nargin < 3 || isempty(tol_rank) || strcmp(tol_rank,'robust')
tol_rank = max(size(X)) * eps(norm(X)); %tolerance level used for rank computations
end
if nargin < 4 || isempty(tol_sv)
tol_sv = 1.e-3; %tolerance level for zero singular value
......
matlab/identification_checks_via_subsets.m
View file @ c4f7c416
......@@ -93,7 +93,11 @@ if ~no_identification_dynamic
if totparam_nbr > modparam_nbr
dDYNAMIC = [zeros(size(ide_dynamic.dDYNAMIC,1),totparam_nbr-modparam_nbr) dDYNAMIC]; %add derivatives wrt stderr and corr parameters
end
rank_dDYNAMIC = rank(full(dDYNAMIC),tol_rank); %compute rank with imposed tolerance level
if strcmp(tol_rank,'robust')
rank_dDYNAMIC = rank(full(dDYNAMIC)); %compute rank with imposed tolerance level
else
rank_dDYNAMIC = rank(full(dDYNAMIC),tol_rank); %compute rank with imposed tolerance level
end
ide_dynamic.rank = rank_dDYNAMIC;
% check rank criteria for full Jacobian
if rank_dDYNAMIC == totparam_nbr
......@@ -112,7 +116,11 @@ end
if ~no_identification_reducedform
dREDUCEDFORM = ide_reducedform.dREDUCEDFORM;
dREDUCEDFORM(ide_reducedform.ind_dREDUCEDFORM,:) = dREDUCEDFORM(ide_reducedform.ind_dREDUCEDFORM,:)./ide_reducedform.norm_dREDUCEDFORM; %normalize
rank_dREDUCEDFORM = rank(full(dREDUCEDFORM),tol_rank); %compute rank with imposed tolerance level
if strcmp(tol_rank,'robust')
rank_dREDUCEDFORM = rank(full(dREDUCEDFORM)); %compute rank with imposed tolerance level
else
rank_dREDUCEDFORM = rank(full(dREDUCEDFORM),tol_rank); %compute rank with imposed tolerance level
end
ide_reducedform.rank = rank_dREDUCEDFORM;
% check rank criteria for full Jacobian
if rank_dREDUCEDFORM == totparam_nbr
......@@ -131,7 +139,11 @@ end
if ~no_identification_moments
dMOMENTS = ide_moments.dMOMENTS;
dMOMENTS(ide_moments.ind_dMOMENTS,:) = dMOMENTS(ide_moments.ind_dMOMENTS,:)./ide_moments.norm_dMOMENTS; %normalize
rank_dMOMENTS = rank(full(dMOMENTS),tol_rank); %compute rank with imposed tolerance level
if strcmp(tol_rank,'robust')
rank_dMOMENTS = rank(full(dMOMENTS)); %compute rank with imposed tolerance level
else
rank_dMOMENTS = rank(full(dMOMENTS),tol_rank); %compute rank with imposed tolerance level
end
ide_moments.rank = rank_dMOMENTS;
% check rank criteria for full Jacobian
if rank_dMOMENTS == totparam_nbr
......@@ -152,7 +164,11 @@ if ~no_identification_spectrum
%alternative normalization
%dSPECTRUM = ide_spectrum.dSPECTRUM;
%dSPECTRUM(ide_spectrum.ind_dSPECTRUM,:) = dSPECTRUM(ide_spectrum.ind_dSPECTRUM,:)./ide_spectrum.norm_dSPECTRUM; %normalize
rank_dSPECTRUM = rank(full(dSPECTRUM),tol_rank); %compute rank with imposed tolerance level
if strcmp(tol_rank,'robust')
rank_dSPECTRUM = rank(full(dSPECTRUM)); %compute rank with imposed tolerance level
else
rank_dSPECTRUM = rank(full(dSPECTRUM),tol_rank); %compute rank with imposed tolerance level
end
ide_spectrum.rank = rank_dSPECTRUM;
% check rank criteria for full Jacobian
if rank_dSPECTRUM == totparam_nbr
......@@ -173,7 +189,11 @@ if ~no_identification_minimal
dMINIMAL(ide_minimal.ind_dMINIMAL,:) = dMINIMAL(ide_minimal.ind_dMINIMAL,:)./ide_minimal.norm_dMINIMAL; %normalize
dMINIMAL_par = dMINIMAL(:,1:totparam_nbr); %part of dMINIMAL that is dependent on parameters
dMINIMAL_rest = dMINIMAL(:,(totparam_nbr+1):end); %part of dMINIMAL that is independent of parameters
rank_dMINIMAL = rank(full(dMINIMAL),tol_rank); %compute rank via SVD see function below
if strcmp(tol_rank,'robust')
rank_dMINIMAL = rank(full(dMINIMAL)); %compute rank via SVD see function below
else
rank_dMINIMAL = rank(full(dMINIMAL),tol_rank); %compute rank via SVD see function below
end
ide_minimal.rank = rank_dMINIMAL;
dMINIMAL_fixed_rank_nbr = size(dMINIMAL_rest,2);
% check rank criteria for full Jacobian
......@@ -193,20 +213,38 @@ end
for j=1:totparam_nbr
if ~no_identification_dynamic
% Columns correspond to single parameters, i.e. full rank would be equal to 1
if rank(full(dDYNAMIC(:,j)),tol_rank) == 0
dynamic_problpars{1} = [dynamic_problpars{1};j];
if strcmp(tol_rank,'robust')
if rank(full(dDYNAMIC(:,j))) == 0
dynamic_problpars{1} = [dynamic_problpars{1};j];
end
else
if rank(full(dDYNAMIC(:,j)),tol_rank) == 0
dynamic_problpars{1} = [dynamic_problpars{1};j];
end
end
end
if ~no_identification_reducedform
% Columns correspond to single parameters, i.e. full rank would be equal to 1
if rank(full(dREDUCEDFORM(:,j)),tol_rank) == 0
reducedform_problpars{1} = [reducedform_problpars{1};j];
if strcmp(tol_rank,'robust')
if rank(full(dREDUCEDFORM(:,j))) == 0
reducedform_problpars{1} = [reducedform_problpars{1};j];
end
else
if rank(full(dREDUCEDFORM(:,j)),tol_rank) == 0
reducedform_problpars{1} = [reducedform_problpars{1};j];
end
end
end
if ~no_identification_moments
if ~no_identification_moments
% Columns correspond to single parameters, i.e. full rank would be equal to 1
if rank(full(dMOMENTS(:,j)),tol_rank) == 0
moments_problpars{1} = [moments_problpars{1};j];
if strcmp(tol_rank,'robust')
if rank(full(dMOMENTS(:,j))) == 0
moments_problpars{1} = [moments_problpars{1};j];
end
else
if rank(full(dMOMENTS(:,j)),tol_rank) == 0
moments_problpars{1} = [moments_problpars{1};j];
end
end
end
if ~no_identification_spectrum
......@@ -218,8 +256,14 @@ for j=1:totparam_nbr
if ~no_identification_minimal
% Columns of dMINIMAL_par correspond to single parameters, needs to be augmented with dMINIMAL_rest (part that is independent of parameters),
% full rank would be equal to 1+dMINIMAL_fixed_rank_nbr
if rank(full([dMINIMAL_par(:,j) dMINIMAL_rest]),tol_rank) == dMINIMAL_fixed_rank_nbr
minimal_problpars{1} = [minimal_problpars{1};j];
if strcmp(tol_rank,'robust')
if rank(full([dMINIMAL_par(:,j) dMINIMAL_rest])) == dMINIMAL_fixed_rank_nbr
minimal_problpars{1} = [minimal_problpars{1};j];
end
else
if rank(full([dMINIMAL_par(:,j) dMINIMAL_rest]),tol_rank) == dMINIMAL_fixed_rank_nbr
minimal_problpars{1} = [minimal_problpars{1};j];
end
end
end
end
......@@ -324,35 +368,55 @@ for j=2:min(length(indtotparam),max_dim_subsets_groups) % Check j-element subset
k_dDYNAMIC = k_dDYNAMIC+1;
if k_dDYNAMIC <= size(indexj_dDYNAMIC,1)
dDYNAMIC_j = dDYNAMIC(:,indexj_dDYNAMIC(k_dDYNAMIC,:)); % pick columns that correspond to parameter subset
rankj_dDYNAMIC(k_dDYNAMIC,1) = rank(full(dDYNAMIC_j),tol_rank); %compute rank with imposed tolerance
if strcmp(tol_rank,'robust')
rankj_dDYNAMIC(k_dDYNAMIC,1) = rank(full(dDYNAMIC_j)); %compute rank with imposed tolerance
else
rankj_dDYNAMIC(k_dDYNAMIC,1) = rank(full(dDYNAMIC_j),tol_rank); %compute rank with imposed tolerance
end
end
end
if ~no_identification_reducedform
k_dREDUCEDFORM = k_dREDUCEDFORM+1;
if k_dREDUCEDFORM <= size(indexj_dREDUCEDFORM,1)
dREDUCEDFORM_j = dREDUCEDFORM(:,indexj_dREDUCEDFORM(k_dREDUCEDFORM,:)); % pick columns that correspond to parameter subset
rankj_dREDUCEDFORM(k_dREDUCEDFORM,1) = rank(full(dREDUCEDFORM_j),tol_rank); %compute rank with imposed tolerance
if strcmp(tol_rank,'robust')
rankj_dREDUCEDFORM(k_dREDUCEDFORM,1) = rank(full(dREDUCEDFORM_j)); %compute rank with imposed tolerance
else
rankj_dREDUCEDFORM(k_dREDUCEDFORM,1) = rank(full(dREDUCEDFORM_j),tol_rank); %compute rank with imposed tolerance
end
end
end
if ~no_identification_moments
k_dMOMENTS = k_dMOMENTS+1;
if k_dMOMENTS <= size(indexj_dMOMENTS,1)
dMOMENTS_j = dMOMENTS(:,indexj_dMOMENTS(k_dMOMENTS,:)); % pick columns that correspond to parameter subset
rankj_dMOMENTS(k_dMOMENTS,1) = rank(full(dMOMENTS_j),tol_rank); %compute rank with imposed tolerance
if strcmp(tol_rank,'robust')
rankj_dMOMENTS(k_dMOMENTS,1) = rank(full(dMOMENTS_j)); %compute rank with imposed tolerance
else
rankj_dMOMENTS(k_dMOMENTS,1) = rank(full(dMOMENTS_j),tol_rank); %compute rank with imposed tolerance
end
end
end
if ~no_identification_minimal
k_dMINIMAL = k_dMINIMAL+1;
if k_dMINIMAL <= size(indexj_dMINIMAL,1)
dMINIMAL_j = [dMINIMAL_par(:,indexj_dMINIMAL(k_dMINIMAL,:)) dMINIMAL_rest]; % pick columns in dMINIMAL_par that correspond to parameter subset and augment with parameter-indepdendet part dMINIMAL_rest
rankj_dMINIMAL(k_dMINIMAL,1) = rank(full(dMINIMAL_j),tol_rank); %compute rank with imposed tolerance
if strcmp(tol_rank,'robust')
rankj_dMINIMAL(k_dMINIMAL,1) = rank(full(dMINIMAL_j)); %compute rank with imposed tolerance
else
rankj_dMINIMAL(k_dMINIMAL,1) = rank(full(dMINIMAL_j),tol_rank); %compute rank with imposed tolerance
end
end
end
if ~no_identification_spectrum
k_dSPECTRUM = k_dSPECTRUM+1;
if k_dSPECTRUM <= size(indexj_dSPECTRUM,1)
dSPECTRUM_j = dSPECTRUM(indexj_dSPECTRUM(k_dSPECTRUM,:),indexj_dSPECTRUM(k_dSPECTRUM,:)); % pick rows and columns that correspond to parameter subset
rankj_dSPECTRUM(k_dSPECTRUM,1) = rank(full(dSPECTRUM_j),tol_rank); % Compute rank with imposed tol
if strcmp(tol_rank,'robust')
rankj_dSPECTRUM(k_dSPECTRUM,1) = rank(full(dSPECTRUM_j)); % Compute rank with imposed tol
else
rankj_dSPECTRUM(k_dSPECTRUM,1) = rank(full(dSPECTRUM_j),tol_rank); % Compute rank with imposed tol
end
end
end
dyn_waitbar(k/maxk,h)
......
matlab/uperm.m
View file @ c4f7c416
% By Willi Mutschler, September 26, 2016. Email: willi@mutschler.eu
% By Bruno Luong
function p = uperm(a)
[u, ~, J] = unique(a);
p = u(up(J, length(a)));
function p = up(J, n)
ktab = histcounts(J,1:max(J));
ktab = histc(J,1:max(J));
l = n;
p = zeros(1, n);
s = 1;
......
tests/identification/as2007/as2007_QT.mod
View file @ c4f7c416
......@@ -2,6 +2,8 @@
% used in their replication file.
% This file is used to check whether the G matrix is computed correctly.
% Created by Willi Mutschler (willi@mutschler.eu)
@#define TOL_RANK = 1e-10
var z g R y pie c piep yp;
varexo e_z e_g e_R ;
parameters tau betta nu phi pistar psi1 psi2 rhor rhog rhoz sig2r sig2g sig2z;
......@@ -72,6 +74,7 @@ identification(parameter_set=calibration,
no_identification_reducedform,
no_identification_moments,
checks_via_subsets=1,
tol_rank=@{TOL_RANK},
max_dim_subsets_groups=4);
load('G_QT'); %note that this is computed using replication files of Qu and Tkachenko (2012)
......@@ -88,7 +91,7 @@ ind_G_QT = (find(max(abs(G_QT'),[],1) > temp.store_options_ident.tol_deriv));
tilda_G_QT = zeros(size(G_QT));
delta_G_QT = sqrt(diag(G_QT(ind_G_QT,ind_G_QT)));
tilda_G_QT(ind_G_QT,ind_G_QT) = G_QT(ind_G_QT,ind_G_QT)./((delta_G_QT)*(delta_G_QT'));
if ~isequal(rank(tilda_G_QT,temp.store_options_ident.tol_rank),rank(tilda_G_dynare,temp.store_options_ident.tol_rank))
if ~isequal(rank(tilda_G_QT,@{TOL_RANK}),rank(tilda_G_dynare,@{TOL_RANK}))
error('ranks are not the same for normalized version')
end
......
tests/identification/as2007/as2007_order_1_2_3.mod
View file @ c4f7c416
......@@ -99,7 +99,7 @@ end;
steady;
check;
stoch_simul(order=3,irf=0,periods=0); %needed for identification(order=3)
stoch_simul(order=3,irf=0); %needed for identification(order=3)
@#for ORDER in [1, 2, 3]
@#for KRONFLAG in [-1, -2, 0]
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment