diff --git a/matlab/disp_identification.m b/matlab/disp_identification.m
index 9c399ebecf5947cc9b39228d18b4be80e2e465b8..716d52f6383ce24231dabffc284f532b28c47c93 100644
--- a/matlab/disp_identification.m
+++ b/matlab/disp_identification.m
@@ -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);
diff --git a/matlab/dynare_identification.m b/matlab/dynare_identification.m
index bb08707df8c379ed1b16278c4c8ddcda9775478f..557083881c239bf1ebd15d776a73958ba23f8c1e 100644
--- a/matlab/dynare_identification.m
+++ b/matlab/dynare_identification.m
@@ -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
diff --git a/matlab/identification_analysis.m b/matlab/identification_analysis.m
index de0be7db44d593cfdba761bff3ac853a225726fb..c03ea8a0ef042f519b06edea6534a3a3be0e0495 100644
--- a/matlab/identification_analysis.m
+++ b/matlab/identification_analysis.m
@@ -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));
diff --git a/matlab/identification_checks.m b/matlab/identification_checks.m
index 5fb971d0fb6d60ab7e340da0befcc083b82b6f4e..cb6717470dd5d1dd87c4efeab975bbdf175da372 100644
--- a/matlab/identification_checks.m
+++ b/matlab/identification_checks.m
@@ -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
diff --git a/matlab/identification_checks_via_subsets.m b/matlab/identification_checks_via_subsets.m
index e3aa35985599def4d1f15e1fe6e8e5c632cb3c27..63a9bbce322bf77d2c7cf30e702f632702237ce1 100644
--- a/matlab/identification_checks_via_subsets.m
+++ b/matlab/identification_checks_via_subsets.m
@@ -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)
diff --git a/matlab/uperm.m b/matlab/uperm.m
index af0d1f6f1f829849bf4e25ca35343f24b83a25e5..bbc690afccd7b511913be92cc270714b72bd5466 100644
--- a/matlab/uperm.m
+++ b/matlab/uperm.m
@@ -1,11 +1,11 @@
-% 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;
diff --git a/tests/identification/as2007/as2007_QT.mod b/tests/identification/as2007/as2007_QT.mod
index 1837bd236601e251159f1d9ee4b6a9b94a605032..715d58db795ffa29294b7ea7848b67af1620ff8f 100644
--- a/tests/identification/as2007/as2007_QT.mod
+++ b/tests/identification/as2007/as2007_QT.mod
@@ -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
diff --git a/tests/identification/as2007/as2007_order_1_2_3.mod b/tests/identification/as2007/as2007_order_1_2_3.mod
index a8fcbfd2d5d4eca3a89de324f277808e8136b5b3..166f36d6939b5af017169de4903e0691f7ab0a47 100644
--- a/tests/identification/as2007/as2007_order_1_2_3.mod
+++ b/tests/identification/as2007/as2007_order_1_2_3.mod
@@ -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]