diff --git a/matlab/dynare_identification.m b/matlab/dynare_identification.m
index 0860935063a3b5600337dfdd2409092056f27f90..b733bb74d9d249b0d26209033498b8829364f5fc 100644
--- a/matlab/dynare_identification.m
+++ b/matlab/dynare_identification.m
@@ -1,4 +1,4 @@
-function [pdraws, STO_REDUCEDFORM, STO_MOMENTS, STO_DYNAMIC, STO_si_dDYNAMIC, STO_si_dREDUCEDFORM, STO_si_dMOMENTS, STO_dSPECTRUM, STO_dMINIMAL] = dynare_identification(options_ident, pdraws0)
+function [pdraws, STO_REDUCEDFORM, STO_MOMENTS, STO_DYNAMIC, STO_si_dDYNAMIC, STO_si_dREDUCEDFORM, STO_si_dMOMENTS, STO_dSPECTRUM, STO_dMINIMAL] = dynare_identification(M_,oo_,options_,bayestopt_,estim_params_,options_ident, pdraws0)
%function [pdraws, STO_REDUCEDFORM, STO_MOMENTS, STO_DYNAMIC, STO_si_dDYNAMIC, STO_si_dREDUCEDFORM, STO_si_dMOMENTS, STO_dSPECTRUM, STO_dMINIMAL] = dynare_identification(options_ident, pdraws0)
% -------------------------------------------------------------------------
% This function is called, when the user specifies identification(...); in the mod file. It prepares all identification analysis:
@@ -65,9 +65,6 @@ function [pdraws, STO_REDUCEDFORM, STO_MOMENTS, STO_DYNAMIC, STO_si_dDYNAMIC, ST
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% =========================================================================
-global M_ options_ oo_ bayestopt_ estim_params_
-
-store_options_ = options_; % store options to restore them at the end
fname = M_.fname; %model name
dname = M_.dname; %model name
@@ -255,11 +252,11 @@ if options_ident.gsa_sample_file
end
pdraws0 = [lpmatx lpmat(istable,:)];
clear lpmat lpmat0 istable;
-elseif nargin==1
+elseif nargin==6
pdraws0=[];
end
external_sample=0;
-if nargin==2 || ~isempty(pdraws0)
+if nargin==7 || ~isempty(pdraws0)
% change settings if there is an external sample provided as input argument
options_ident.prior_mc = size(pdraws0,1);
options_ident.load_ident_files = 0;
@@ -514,8 +511,8 @@ if iload <=0
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 && ~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, M_.fname, options_, estim_params_, parameters_TeX, name_tex);
+ plot_identification(M_,params, ide_moments_point, ide_hess_point, ide_reducedform_point, ide_dynamic_point, options_ident.advanced, parameters, name, ...
+ IdentifDirectoryName, M_.fname, options_, estim_params_, bayestopt_, parameters_TeX, name_tex);
end
if SampleSize > 1
@@ -864,8 +861,8 @@ if iload
disp_identification(ide_hess_point.params, ide_reducedform_point, ide_moments_point, ide_spectrum_point, ide_minimal_point, name, options_ident);
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, M_.fname, options_, estim_params_, [], name_tex);
+ plot_identification(M_,ide_hess_point.params, ide_moments_point, ide_hess_point, ide_reducedform_point, ide_dynamic_point, options_ident.advanced, parameters, name, ...
+ IdentifDirectoryName, M_.fname, options_, estim_params_, bayestopt_, [], name_tex);
end
end
@@ -879,8 +876,8 @@ if SampleSize > 1
options_ident.advanced = advanced0; % reset advanced setting
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, M_.fname, options_, estim_params_, [], name_tex);
+ plot_identification(M_, pdraws, IDE_MOMENTS, ide_hess_point, IDE_REDUCEDFORM, IDE_DYNAMIC, options_ident.advanced, 'MC sample ', name, ...
+ IdentifDirectoryName, M_.fname, options_, estim_params_, bayestopt_, [], name_tex);
end
%advanced display and plots for MC Sample, i.e. look at draws with highest/lowest condition number
if options_ident.advanced
@@ -905,8 +902,8 @@ if SampleSize > 1
options_ident.advanced = advanced0; %reset advanced setting
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, M_.fname, options_, estim_params_, tittxt, name_tex);
+ plot_identification(M_, pdraws(jmax,:), ide_moments_max, ide_hess_max, ide_reducedform_max, ide_dynamic_max, 1, tittxt, name, ...
+ IdentifDirectoryName, M_.fname, options_, estim_params_, bayestopt_, tittxt, name_tex);
end
% SMALLEST condition number
@@ -924,8 +921,8 @@ if SampleSize > 1
options_ident.advanced = advanced0; %reset advanced setting
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, M_.fname, options_, estim_params_, tittxt,name_tex);
+ plot_identification(M_, pdraws(jmin,:),ide_moments_min,ide_hess_min,ide_reducedform_min,ide_dynamic_min,1,tittxt,name,...
+ IdentifDirectoryName, M_.fname, options_, estim_params_, bayestopt_, tittxt,name_tex);
end
% reset nodisplay option
options_.nodisplay = store_nodisplay;
@@ -946,8 +943,8 @@ if SampleSize > 1
options_ident.advanced = advanced0; % reset advanced
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, M_.fname, options_, estim_params_, tittxt, name_tex);
+ plot_identification(M_, pdraws(jcrit(j),:), ide_moments_(j), ide_hess_(j), ide_reducedform_(j), ide_dynamic_(j), 1, tittxt, name, ...
+ IdentifDirectoryName, M_.fname, options_, estim_params_, bayestopt_, tittxt, name_tex);
end
end
if ~iload
@@ -963,6 +960,4 @@ end
%reset warning state
warning_config;
-fprintf('\n==== Identification analysis completed ====\n\n')
-
-options_ = store_options_; %restore options set
+fprintf('\n==== Identification analysis completed ====\n\n')
\ No newline at end of file
diff --git a/matlab/dynare_sensitivity.m b/matlab/dynare_sensitivity.m
index 1975def14205ad120b2557566394b7a3f48be1f3..405fb92f1d517c0053561da4e0c8d76a4baff0c6 100644
--- a/matlab/dynare_sensitivity.m
+++ b/matlab/dynare_sensitivity.m
@@ -1,10 +1,17 @@
-function x0=dynare_sensitivity(options_gsa)
+function x0=dynare_sensitivity(M_,oo_,options_,bayestopt_,estim_params_,options_gsa)
% Frontend to the Sensitivity Analysis Toolbox for DYNARE
+% Inputs:
+% - M_ [structure] Matlab's structure describing the model
+% - oo_ [structure] Matlab's structure describing the results
+% - options_ [structure] Matlab's structure describing the current options
+% - bayestopt_ [structure] describing the priors
+% - estim_params_ [structure] characterizing parameters to be estimated
+% - options_gsa [structure] Matlab's structure describing the GSA options
%
% Reference:
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
-% Copyright © 2008-2018 Dynare Team
+% Copyright © 2008-2023 Dynare Team
%
% This file is part of Dynare.
%
@@ -21,8 +28,6 @@ function x0=dynare_sensitivity(options_gsa)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-global M_ options_ oo_ bayestopt_ estim_params_
-
if options_.dsge_var
error('Identification does not support DSGE-VARs at the current stage')
end
@@ -87,9 +92,6 @@ if options_.order~=1
options_.order = 1;
end
-original_prior_trunc = options_.prior_trunc;
-original_qz_criterium = options_.qz_criterium;
-
if ~isempty(options_gsa.datafile) || isempty(bayestopt_) || options_gsa.rmse
if isempty(options_gsa.datafile) && options_gsa.rmse
disp('The data file and all relevant estimation options ')
@@ -208,9 +210,9 @@ options_gsa = set_default_option(options_gsa,'var_rmse', options_.varobs);
options_gsa.var_rmse_tex={};
for ii=1:length(options_gsa.var_rmse)
temp_name = M_.endo_names_tex{strmatch(options_gsa.var_rmse{ii}, M_.endo_names, 'exact')};
- options_gsa.var_rmse_tex = vertcat(options_gsa.var_rmse_tex, temp_name);
+ options_gsa.var_rmse_tex = vertcat(options_gsa.var_rmse_tex, ['$' temp_name '$']);
end
-options_gsa.var_rmse_tex = cellfun(@(x) horzcat('$', x, '$'), options_gsa.var_rmse_tex, 'UniformOutput', false);
+options_gsa.varobs_tex = cellfun(@(x) horzcat('$', x, '$'), M_.endo_names_tex(options_.varobs_id), 'UniformOutput', false);
options_gsa = set_default_option(options_gsa,'pfilt_rmse', 0.1);
options_gsa = set_default_option(options_gsa,'istart_rmse', options_.presample+1);
options_gsa = set_default_option(options_gsa,'alpha_rmse', 0.001);
@@ -281,7 +283,7 @@ if (options_gsa.load_stab || options_gsa.load_rmse || options_gsa.load_redform)
end
if options_gsa.stab && ~options_gsa.ppost
- x0 = stab_map_(OutputDirectoryName,options_gsa);
+ x0 = stab_map_(OutputDirectoryName,options_gsa,M_,oo_,options_,bayestopt_,estim_params_);
if isempty(x0)
skipline()
disp('Sensitivity computations stopped: no parameter set provided a unique solution')
@@ -295,7 +297,7 @@ if ~isempty(options_gsa.moment_calibration) || ~isempty(options_gsa.irf_calibrat
end
if options_gsa.identification
- map_ident_(OutputDirectoryName,options_gsa,M_,options_,bayestopt_,oo_.dr,estim_params_.param_vals);
+ map_ident_(OutputDirectoryName,options_gsa,M_,oo_,options_,estim_params_,bayestopt_);
end
if options_gsa.redform && ~isempty(options_gsa.namendo)
@@ -303,7 +305,7 @@ if options_gsa.redform && ~isempty(options_gsa.namendo)
filnam = dir([M_.dname filesep 'metropolis' filesep '*param_irf*.mat']);
lpmat=[];
for j=1:length(filnam)
- load ([M_.dname filesep 'metropolis' filesep M_.fname '_param_irf' int2str(j) '.mat'])
+ load ([M_.dname filesep 'metropolis' filesep M_.fname '_param_irf' int2str(j) '.mat'],'stock')
lpmat=[lpmat; stock];
end
clear stock
@@ -321,7 +323,7 @@ if options_gsa.redform && ~isempty(options_gsa.namendo)
save([OutputDirectoryName filesep M_.fname '_mc.mat'],'lpmat','lpmat0','istable','iunstable','iwrong','iindeterm')
options_gsa.load_stab=1;
- x0 = stab_map_(OutputDirectoryName,options_gsa);
+ x0 = stab_map_(OutputDirectoryName,options_gsa,M_,oo_,options_,bayestopt_,estim_params_);
end
if strmatch(':',options_gsa.namendo,'exact')
options_gsa.namendo = M_.endo_names(1:M_.orig_endo_nbr);
@@ -332,9 +334,8 @@ if options_gsa.redform && ~isempty(options_gsa.namendo)
if strmatch(':',options_gsa.namlagendo,'exact')
options_gsa.namlagendo = M_.endo_names(1:M_.orig_endo_nbr);
end
- % options_.opt_gsa = options_gsa;
if options_gsa.morris==1
- redform_screen(OutputDirectoryName,options_gsa, estim_params_.param_vals, M_, oo_.dr, options_, bayestopt_);
+ redform_screen(OutputDirectoryName,options_gsa, estim_params_, M_, oo_.dr, options_, bayestopt_);
else
% check existence of the SS_ANOVA toolbox
if isempty(options_gsa.threshold_redform) && ~(exist('gsa_sdp','file')==6 || exist('gsa_sdp','file')==2)
@@ -345,7 +346,7 @@ if options_gsa.redform && ~isempty(options_gsa.namendo)
fprintf('After obtaining the files, you need to unpack them and set a Matlab Path to those files.\n')
error('SS-ANOVA-R Toolbox missing!')
end
- redform_map(OutputDirectoryName,options_gsa,M_,estim_params_,options_,bayestopt_,oo_.dr);
+ redform_map(OutputDirectoryName,options_gsa,M_,estim_params_,options_,bayestopt_,oo_);
end
end
% RMSE mapping
@@ -403,9 +404,6 @@ if options_gsa.rmse
filt_mc_(OutputDirectoryName,options_gsa,dataset_,dataset_info,M_,oo_,options_,bayestopt_,estim_params_);
end
options_.opt_gsa = options_gsa;
-options_.prior_trunc=original_prior_trunc;
-options_.qz_criterium=original_qz_criterium ;
-
if options_gsa.glue
dr_ = oo_.dr;
@@ -444,7 +442,6 @@ if options_gsa.glue
vj = options_.varobs{j};
jxj = strmatch(vj,M_.endo_names(dr_.order_var),'exact');
- js = strmatch(vj,M_.endo_names,'exact');
if ~options_gsa.ppost
Out(j).data = jxj;
Out(j).time = [pwd,'/',OutputDirectoryName];
diff --git a/matlab/gsa/filt_mc_.m b/matlab/gsa/filt_mc_.m
index 06c19404d4c827143f271253490b06761f641c1f..1d7e0944315bd893730ca0700eae189cce9ca044 100644
--- a/matlab/gsa/filt_mc_.m
+++ b/matlab/gsa/filt_mc_.m
@@ -84,11 +84,15 @@ if ~options_.nograph
end
for j=1:length(a)
if strmatch([fname_,tmp],a(j).name)
- disp(a(j).name)
+ if options_.debug
+ disp(a(j).name)
+ end
delete([OutDir,filesep,a(j).name])
end
if strmatch([fname_,tmp1],a(j).name)
- disp(a(j).name)
+ if options_.debug
+ disp(a(j).name)
+ end
delete([OutDir,filesep,a(j).name])
end
end
@@ -306,7 +310,7 @@ if ~options_.opt_gsa.ppost && options_.opt_gsa.lik_only
options_mcf.title = atitle;
options_mcf.beha_title = 'better posterior kernel';
options_mcf.nobeha_title = 'worse posterior kernel';
- mcf_analysis(x, ipost(1:nfilt), ipost(nfilt+1:end), options_mcf, options_, bayestopt_, estim_params_);
+ mcf_analysis(x, ipost(1:nfilt), ipost(nfilt+1:end), options_mcf, M_, options_, bayestopt_, estim_params_);
if options_.opt_gsa.pprior
anam = 'rmse_prior_lik';
atitle = 'RMSE prior: Log Likelihood Kernel';
@@ -319,7 +323,7 @@ if ~options_.opt_gsa.ppost && options_.opt_gsa.lik_only
options_mcf.title = atitle;
options_mcf.beha_title = 'better likelihood';
options_mcf.nobeha_title = 'worse likelihood';
- mcf_analysis(x, ilik(1:nfilt), ilik(nfilt+1:end), options_mcf, options_, bayestopt_, estim_params_);
+ mcf_analysis(x, ilik(1:nfilt), ilik(nfilt+1:end), options_mcf, M_, options_, bayestopt_, estim_params_);
else
if options_.opt_gsa.ppost
@@ -729,10 +733,15 @@ else
for iy = 1:length(vvarvecm)
options_mcf.amcf_name = [asname '_' vvarvecm{iy} '_map' ];
options_mcf.amcf_title = [atitle ' ' vvarvecm{iy}];
- options_mcf.beha_title = ['better fit of ' vvarvecm{iy}];
- options_mcf.nobeha_title = ['worse fit of ' vvarvecm{iy}];
+ if options_.TeX
+ options_mcf.beha_title = ['better fit of ' vvarvecm_tex{iy}];
+ options_mcf.nobeha_title = ['worse fit of ' vvarvecm_tex{iy}];
+ else
+ options_mcf.beha_title = ['better fit of ' vvarvecm{iy}];
+ options_mcf.nobeha_title = ['worse fit of ' vvarvecm{iy}];
+ end
options_mcf.title = ['the fit of ' vvarvecm{iy}];
- mcf_analysis(x, ixx(1:nfilt0(iy),iy), ixx(nfilt0(iy)+1:end,iy), options_mcf, options_, bayestopt_, estim_params_);
+ mcf_analysis(x, ixx(1:nfilt0(iy),iy), ixx(nfilt0(iy)+1:end,iy), options_mcf, M_, options_, bayestopt_, estim_params_);
end
for iy = 1:length(vvarvecm)
ipar = find(any(squeeze(PPV(iy,:,:))<alpha));
@@ -877,7 +886,7 @@ pnames_tex=cell(np,1);
for ii=1:length(bayestopt_.name)
if options_.TeX
[param_name_temp, param_name_tex_temp]= get_the_name(ii,options_.TeX,M_,estim_params_,options_);
- pnames_tex{ii,1} = strrep(param_name_tex_temp,'$','');
+ pnames_tex{ii,1} = param_name_tex_temp;
pnames{ii,1} = param_name_temp;
else
param_name_temp = get_the_name(ii,options_.TeX,M_,estim_params_,options_);
diff --git a/matlab/gsa/map_calibration.m b/matlab/gsa/map_calibration.m
index 0fce0c6fdab6f5fa2d2d97e8237defc60dc396bd..395544e00a03a8bced09d520c2df7f2bf5f9b6a8 100644
--- a/matlab/gsa/map_calibration.m
+++ b/matlab/gsa/map_calibration.m
@@ -253,7 +253,7 @@ if ~isempty(indx_irf)
options_mcf.nobeha_title = 'NO IRF restriction';
options_mcf.title = atitle0;
if ~isempty(indx1) && ~isempty(indx2)
- mcf_analysis(xmat(:,nshock+1:end), indx1, indx2, options_mcf, options_, bayestopt_, estim_params_);
+ mcf_analysis(xmat(:,nshock+1:end), indx1, indx2, options_mcf, M_, options_, bayestopt_, estim_params_);
end
end
for ij=1:nbr_irf_couples
@@ -305,7 +305,7 @@ if ~isempty(indx_irf)
options_mcf.nobeha_title = 'NO IRF restriction';
options_mcf.title = atitle0;
if ~isempty(indx1) && ~isempty(indx2)
- mcf_analysis(xmat(:,nshock+1:end), indx1, indx2, options_mcf, options_, bayestopt_, estim_params_);
+ mcf_analysis(xmat(:,nshock+1:end), indx1, indx2, options_mcf, M_, options_, bayestopt_, estim_params_);
end
end
end
@@ -354,7 +354,7 @@ if ~isempty(indx_moment)
for jj=1:np
if options_.TeX
[param_name_temp, param_name_tex_temp]= get_the_name(jj,options_.TeX,M_,estim_params_,options_);
- name_tex{jj,1} = strrep(param_name_tex_temp,'$','');
+ name_tex{jj,1} = param_name_tex_temp;
name{jj,1} = param_name_temp;
else
param_name_temp = get_the_name(jj,options_.TeX,M_,estim_params_,options_);
@@ -448,7 +448,7 @@ if ~isempty(indx_moment)
options_mcf.nobeha_title = 'NO moment restriction';
options_mcf.title = atitle0;
if ~isempty(indx1) && ~isempty(indx2)
- mcf_analysis(xmat, indx1, indx2, options_mcf, options_, bayestopt_, estim_params_);
+ mcf_analysis(xmat, indx1, indx2, options_mcf, M_, options_, bayestopt_, estim_params_);
end
end
for ij=1:nbr_moment_couples
@@ -501,7 +501,7 @@ if ~isempty(indx_moment)
options_mcf.nobeha_title = 'NO moment restriction';
options_mcf.title = atitle0;
if ~isempty(indx1) && ~isempty(indx2)
- mcf_analysis(xmat, indx1, indx2, options_mcf, options_, bayestopt_, estim_params_);
+ mcf_analysis(xmat, indx1, indx2, options_mcf, M_, options_, bayestopt_, estim_params_);
end
end
end
diff --git a/matlab/gsa/map_ident_.m b/matlab/gsa/map_ident_.m
index cb46f4ff95897e505fa29f4e3869f939061e0cff..cabe148f33ca09a5c7242de72985f02622fef1c7 100644
--- a/matlab/gsa/map_ident_.m
+++ b/matlab/gsa/map_ident_.m
@@ -1,13 +1,13 @@
-function map_ident_(OutputDirectoryName,opt_gsa,M_,options_,bayestopt_,dr,param_vals)
-% map_ident_(OutputDirectoryName,opt_gsa,M_,oo_,options_,bayestopt_,param_vals)
+function map_ident_(OutputDirectoryName,opt_gsa,M_,oo_,options_,estim_params_,bayestopt_)
+% map_ident_(OutputDirectoryName,opt_gsa,M_,oo_,options_,estim_params_,bayestopt_)
% Inputs
% - OutputDirectoryName [string] name of the output directory
% - opt_gsa [structure] GSA options structure
% - M_ [structure] Matlab's structure describing the model
+% - oo_ [structure] Matlab's structure describing the results
% - options_ [structure] Matlab's structure describing the current options
+% - estim_params_ [structure] characterizing parameters to be estimated
% - bayestopt_ [structure] describing the priors
-% - dr [structure] decision rules
-% - param_vals [structure] characterizing parameters to be estimated
% Written by Marco Ratto
% Joint Research Centre, The European Commission,
@@ -32,12 +32,13 @@ function map_ident_(OutputDirectoryName,opt_gsa,M_,options_,bayestopt_,dr,param_
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
fname_ = M_.fname;
+dr=oo_.dr;
nliv = opt_gsa.morris_nliv;
itrans = opt_gsa.trans_ident;
-np = size(param_vals,1);
+np = size(estim_params_.param_vals,1);
-pnames = M_.param_names(param_vals(:,1));
+pnames = M_.param_names(estim_params_.param_vals(:,1));
if opt_gsa.pprior
filetoload=[OutputDirectoryName '/' fname_ '_prior'];
else
@@ -60,15 +61,15 @@ if opt_gsa.load_ident_files==0
mss = teff(mss(:,istable),Nsam,istable);
yys = teff(yys(dr.order_var,istable),Nsam,istable);
if exist('T','var')
- [vdec, cc, ac] = mc_moments(T, lpmatx, dr, M_, options_);
+ [vdec, cc, ac] = mc_moments(T, lpmatx, dr, M_, options_, estim_params_);
else
return
end
if opt_gsa.morris==2
- pdraws = dynare_identification(options_.options_ident,[lpmatx lpmat(istable,:)]);
+ pdraws = dynare_identification(M_,oo_,options_,bayestopt_,estim_params_,options_.options_ident,[lpmatx lpmat(istable,:)]);
if ~isempty(pdraws) && max(max(abs(pdraws-[lpmatx lpmat(istable,:)])))==0
- disp(['Sample check OK ', num2str(max(max(abs(pdraws-[lpmatx lpmat(istable,:)]))))]),
+ disp(['Sample check OK. Largest difference: ', num2str(max(max(abs(pdraws-[lpmatx lpmat(istable,:)]))))]),
clear pdraws;
end
clear GAM gas
@@ -88,7 +89,11 @@ if opt_gsa.load_ident_files==0
set(gca,'xlim',[0.5 size(options_.varobs,1)+0.5])
set(gca,'ylim',[-2 102])
for ip=1:size(options_.varobs,1)
- text(ip,-4,deblank(options_.varobs(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
+ if options_.TeX
+ text(ip,-4,deblank(opt_gsa.varobs_tex(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','latex')
+ else
+ text(ip,-4,deblank(options_.varobs(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
+ end
end
xlabel(' ')
ylabel(' ')
@@ -118,7 +123,7 @@ if opt_gsa.load_ident_files==0
[Aa,Bb] = kalman_transition_matrix(dr,iv,ic);
A = zeros(size(Aa,1),size(Aa,2)+size(Aa,1),length(istable));
if ~isempty(lpmatx)
- set_shocks_param(lpmatx(1,:));
+ M_=set_shocks_param(M_,estim_params_,lpmatx(1,:));
end
A(:,:,1)=[Aa, triu(Bb*M_.Sigma_e*Bb')];
for j=2:length(istable)
@@ -126,7 +131,7 @@ if opt_gsa.load_ident_files==0
dr.ghu = T(:, (nc1-M_.exo_nbr+1):end, j);
[Aa,Bb] = kalman_transition_matrix(dr, iv, ic);
if ~isempty(lpmatx)
- set_shocks_param(lpmatx(j,:));
+ M_=set_shocks_param(M_,estim_params_,lpmatx(j,:));
end
A(:,:,j)=[Aa, triu(Bb*M_.Sigma_e*Bb')];
end
@@ -166,7 +171,12 @@ if opt_gsa.morris==1
set(gca,'ylim',[0 ydum(2)])
set(gca,'position',[0.13 0.2 0.775 0.7])
for ip=1:npT
- text(ip,-2,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
+ if options_.TeX
+ [~, param_name_tex_temp]= get_the_name(ip,options_.TeX,M_,estim_params_,options_);
+ text(ip,-2,param_name_tex_temp,'rotation',90,'HorizontalAlignment','right','interpreter','latex')
+ else
+ text(ip,-2,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
+ end
end
xlabel(' ')
title('Elementary effects variance decomposition')
@@ -198,7 +208,12 @@ if opt_gsa.morris==1
set(gca,'ylim',[0 1])
set(gca,'position',[0.13 0.2 0.775 0.7])
for ip=1:npT
- text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
+ if options_.TeX
+ [~, param_name_tex_temp]= get_the_name(ip,options_.TeX,M_,estim_params_,options_);
+ text(ip,-0.02,param_name_tex_temp,'rotation',90,'HorizontalAlignment','right','interpreter','latex')
+ else
+ text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
+ end
end
xlabel(' ')
title('Elementary effects in the moments')
@@ -241,7 +256,12 @@ if opt_gsa.morris==1
set(gca,'position',[0.13 0.2 0.775 0.7])
xlabel(' ')
for ip=1:npT
- text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
+ if options_.TeX
+ [~, param_name_tex_temp]= get_the_name(ip,options_.TeX,M_,estim_params_,options_);
+ text(ip,-0.02,param_name_tex_temp,'rotation',90,'HorizontalAlignment','right','interpreter','latex')
+ else
+ text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
+ end
end
xlabel(' ')
title('Elementary effects in the model')
diff --git a/matlab/gsa/mc_moments.m b/matlab/gsa/mc_moments.m
index b942368bc500b3c10cdc661a51cec9e611f1bb94..d32d01d1316a09a7e6dd6e690302bc6eccac925a 100644
--- a/matlab/gsa/mc_moments.m
+++ b/matlab/gsa/mc_moments.m
@@ -1,10 +1,11 @@
-function [vdec, cc, ac] = mc_moments(mm, ss, dr, M_, options_)
-% [vdec, cc, ac] = mc_moments(mm, ss, dr, M_, options_)
+function [vdec, cc, ac] = mc_moments(mm, ss, dr, M_, options_, estim_params_)
+% [vdec, cc, ac] = mc_moments(mm, ss, dr, M_, options_,estim_params_)
% Conduct Monte Carlo simulation of second moments for GSA
% Inputs:
% - dr [structure] decision rules
% - M_ [structure] model structure
% - options_ [structure] Matlab's structure describing the current options
+% - estim_params_ [structure] characterizing parameters to be estimated
%
% Outputs:
% - vdec [double] variance decomposition matrix
@@ -41,7 +42,7 @@ for j=1:nsam
dr.ghx = mm(:, 1:(nc1-M_.exo_nbr),j);
dr.ghu = mm(:, (nc1-M_.exo_nbr+1):end, j);
if ~isempty(ss)
- set_shocks_param(ss(j,:));
+ M_=set_shocks_param(M_,estim_params_,ss(j,:));
end
[vdec(:,:,j), corr, autocorr] = th_moments(dr,options_,M_);
cc(:,:,j)=triu(corr);
diff --git a/matlab/gsa/mcf_analysis.m b/matlab/gsa/mcf_analysis.m
index c901ef2c36ed3e0668e38104208f317ab9369d26..916013721ef224417b79c89c21df3ed14affe65f 100644
--- a/matlab/gsa/mcf_analysis.m
+++ b/matlab/gsa/mcf_analysis.m
@@ -1,10 +1,11 @@
-function indmcf = mcf_analysis(lpmat, ibeha, inobeha, options_mcf, options_, bayestopt_, estim_params_)
-% indmcf = mcf_analysis(lpmat, ibeha, inobeha, options_mcf, options_, bayestopt_, estim_params_)
+function indmcf = mcf_analysis(lpmat, ibeha, inobeha, options_mcf, M_, options_, bayestopt_, estim_params_)
+% indmcf = mcf_analysis(lpmat, ibeha, inobeha, options_mcf, M_, options_, bayestopt_, estim_params_)
% Inputs:
% - lpmat [double] Monte Carlo matrix
% - ibeha [integer] index of behavioural runs
% - inobeha [integer] index of non-behavioural runs
% - options_gsa_ [structure] GSA options_
+% - M_ [structure] describing the model
% - options_ [structure] describing the options
% - bayestopt_ [structure] describing the priors
% - estim_params_ [structure] characterizing parameters to be estimated
@@ -81,8 +82,8 @@ if ~isempty(indmcf)
end
if length(ibeha)>10 && length(inobeha)>10
- indcorr1 = stab_map_2(lpmat(ibeha,:),alpha2, pvalue_corr, beha_title);
- indcorr2 = stab_map_2(lpmat(inobeha,:),alpha2, pvalue_corr, nobeha_title);
+ indcorr1 = stab_map_2(lpmat(ibeha,:),alpha2, pvalue_corr, M_, options_, bayestopt_, estim_params_, beha_title);
+ indcorr2 = stab_map_2(lpmat(inobeha,:),alpha2, pvalue_corr, M_, options_, bayestopt_, estim_params_, nobeha_title);
indcorr = union(indcorr1(:), indcorr2(:));
indcorr = indcorr(~ismember(indcorr(:),indmcf));
indmcf = [indmcf(:); indcorr(:)];
@@ -90,7 +91,9 @@ end
if ~isempty(indmcf) && ~options_.nograph
skipline()
xx=[];
- if ~ isempty(xparam1), xx=xparam1(indmcf); end
+ if ~ isempty(xparam1)
+ xx=xparam1(indmcf);
+ end
scatter_mcf(lpmat(ibeha,indmcf),lpmat(inobeha,indmcf), param_names_tex(indmcf), ...
'.', [fname_,'_',amcf_name], OutputDirectoryName, amcf_title,xx, options_, ...
beha_title, nobeha_title)
diff --git a/matlab/gsa/prior_draw_gsa.m b/matlab/gsa/prior_draw_gsa.m
index 85a26f3ca8cd2fa9408bf2d14f9b2a1f595ae545..58731ec0a1120d248ab356e94a8f1ee464f4136d 100644
--- a/matlab/gsa/prior_draw_gsa.m
+++ b/matlab/gsa/prior_draw_gsa.m
@@ -1,9 +1,13 @@
-function pdraw = prior_draw_gsa(init,rdraw)
+function pdraw = prior_draw_gsa(M_,bayestopt_,options_,estim_params_,init,rdraw)
% Draws from the prior distributions for use with Sensitivity Toolbox for DYNARE
%
% INPUTS
-% o init [integer] scalar equal to 1 (first call) or 0.
-% o rdraw
+% - M_ [structure] describing the model
+% - bayestopt_ [structure] describing the priors
+% - options_ [structure] describing the options
+% - estim_params_ [structure] characterizing parameters to be estimated
+% - init [integer] scalar equal to 1 (first call) or 0.
+% - rdraw
%
% OUTPUTS
% o pdraw [double] draw from the joint prior density.
@@ -35,8 +39,7 @@ function pdraw = prior_draw_gsa(init,rdraw)
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-
-global bayestopt_ options_ estim_params_ M_
+
persistent npar pshape p6 p7 p3 p4 lbcum ubcum
if init
@@ -49,7 +52,7 @@ if init
pdraw = zeros(npar,1);
lbcum = zeros(npar,1);
ubcum = ones(npar,1);
- [~,~,~,lb,ub,~] = set_prior(estim_params_,M_,options_); %Prepare bounds
+ [~,~,~,lb,ub] = set_prior(estim_params_,M_,options_); %Prepare bounds
if ~isempty(bayestopt_) && any(bayestopt_.pshape > 0)
% Set prior bounds
bounds = prior_bounds(bayestopt_, options_.prior_trunc);
@@ -64,29 +67,29 @@ if init
% set bounds for cumulative probabilities
for i = 1:npar
switch pshape(i)
- case 5% Uniform prior.
- p4(i) = min(p4(i),bounds.ub(i));
- p3(i) = max(p3(i),bounds.lb(i));
- case 3% Gaussian prior.
- lbcum(i) = 0.5 * erfc(-(bounds.lb(i)-p6(i))/p7(i) ./ sqrt(2));
- ubcum(i) = 0.5 * erfc(-(bounds.ub(i)-p6(i))/p7(i) ./ sqrt(2));
- case 2% Gamma prior.
- lbcum(i) = gamcdf(bounds.lb(i)-p3(i),p6(i),p7(i));
- ubcum(i) = gamcdf(bounds.ub(i)-p3(i),p6(i),p7(i));
case 1% Beta distribution (TODO: generalized beta distribution)
lbcum(i) = betainc((bounds.lb(i)-p3(i))./(p4(i)-p3(i)),p6(i),p7(i));
ubcum(i) = betainc((bounds.ub(i)-p3(i))./(p4(i)-p3(i)),p6(i),p7(i));
+ case 2% Gamma prior.
+ lbcum(i) = gamcdf(bounds.lb(i)-p3(i),p6(i),p7(i));
+ ubcum(i) = gamcdf(bounds.ub(i)-p3(i),p6(i),p7(i));
+ case 3% Gaussian prior.
+ lbcum(i) = 0.5 * erfc(-(bounds.lb(i)-p6(i))/p7(i) ./ sqrt(2));
+ ubcum(i) = 0.5 * erfc(-(bounds.ub(i)-p6(i))/p7(i) ./ sqrt(2));
case 4% INV-GAMMA1 distribution
% TO BE CHECKED
lbcum(i) = gamcdf(1/(bounds.ub(i)-p3(i))^2,p7(i)/2,2/p6(i));
ubcum(i) = gamcdf(1/(bounds.lb(i)-p3(i))^2,p7(i)/2,2/p6(i));
+ case 5% Uniform prior.
+ p4(i) = min(p4(i),bounds.ub(i));
+ p3(i) = max(p3(i),bounds.lb(i));
case 6% INV-GAMMA2 distribution
% TO BE CHECKED
lbcum(i) = gamcdf(1/(bounds.ub(i)-p3(i)),p7(i)/2,2/p6(i));
ubcum(i) = gamcdf(1/(bounds.lb(i)-p3(i)),p7(i)/2,2/p6(i));
case 8
- lbcum(i) = weibcdf(bounds.lb(i)-p3(i),p6(i),p7(i));
- ubcum(i) = weibcdf(bounds.ub(i)-p3(i),p6(i),p7(i));
+ lbcum(i) = wblcdf(bounds.lb(i)-p3(i),p6(i),p7(i));
+ ubcum(i) = wblcdf(bounds.ub(i)-p3(i),p6(i),p7(i));
otherwise
% Nothing to do here.
end
@@ -94,7 +97,7 @@ if init
return
end
-
+pdraw=NaN(size(rdraw,1),npar);
for i = 1:npar
rdraw(:,i) = rdraw(:,i).*(ubcum(i)-lbcum(i))+lbcum(i);
switch pshape(i)
diff --git a/matlab/gsa/redform_map.m b/matlab/gsa/redform_map.m
index 221a052c091a09ae89d57eb43dba1a3fe77c28b9..d8681ce2da2cf9c5caa7ec70d0e6eda6cff4e96c 100644
--- a/matlab/gsa/redform_map.m
+++ b/matlab/gsa/redform_map.m
@@ -1,5 +1,5 @@
-function redform_map(dirname,options_gsa_,M_,estim_params_,options_,bayestopt_,dr)
-% redform_map(dirname,options_gsa_,M_,estim_params_,options_,bayestopt_,dr)
+function redform_map(dirname,options_gsa_,M_,estim_params_,options_,bayestopt_,oo_)
+% redform_map(dirname,options_gsa_,M_,estim_params_,options_,bayestopt_,oo_)
% Inputs:
% - dirname [string] name of the output directory
% - options_gsa_ [structure] GSA options_
@@ -7,7 +7,7 @@ function redform_map(dirname,options_gsa_,M_,estim_params_,options_,bayestopt_,d
% - estim_params_ [structure] characterizing parameters to be estimated
% - options_ [structure] describing the options
% - bayestopt_ [structure] describing the priors
-% - dr [structure] storing the decision rules
+% - oo_ [structure] storing the results
%
% Written by Marco Ratto
% Joint Research Centre, The European Commission,
@@ -82,7 +82,7 @@ options_mcf.fname_ = M_.fname;
options_mcf.OutputDirectoryName = adir;
if ~exist('T','var')
- stab_map_(dirname,options_gsa_);
+ stab_map_(dirname,options_gsa_,M_,oo_,options_,bayestopt_,estim_params_);
if pprior
load([dirname,filesep,M_.fname,'_prior'],'T');
else
@@ -134,7 +134,7 @@ lpmat0=[];
js=0;
for j = 1:length(anamendo)
namendo = anamendo{j};
- iendo = strmatch(namendo, M_.endo_names(dr.order_var), 'exact');
+ iendo = strmatch(namendo, M_.endo_names(oo_.dr.order_var), 'exact');
ifig = 0;
iplo = 0;
for jx = 1:length(anamexo)
@@ -205,7 +205,7 @@ for j = 1:length(anamendo)
options_mcf.OutputDirectoryName = xdir;
if ~isempty(iy) && ~isempty(iyc)
fprintf(['%4.1f%% of the ',type,' support matches ',atitle0,'\n'],length(iy)/length(y0)*100)
- icheck = mcf_analysis(x0, iy, iyc, options_mcf, options_, bayestopt_, estim_params_);
+ icheck = mcf_analysis(x0, iy, iyc, options_mcf, M_, options_, bayestopt_, estim_params_);
lpmat=x0(iy,:);
if nshocks
@@ -289,7 +289,7 @@ for j = 1:length(anamendo)
iplo=0;
for je=1:length(anamlagendo)
namlagendo = anamlagendo{je};
- ilagendo=strmatch(namlagendo, M_.endo_names(dr.order_var(M_.nstatic+1:M_.nstatic+nsok)), 'exact');
+ ilagendo=strmatch(namlagendo, M_.endo_names(oo_.dr.order_var(M_.nstatic+1:M_.nstatic+nsok)), 'exact');
skipline()
disp(['[', namendo,' vs lagged ',namlagendo,']'])
@@ -356,7 +356,7 @@ for j = 1:length(anamendo)
if ~isempty(iy) && ~isempty(iyc)
fprintf(['%4.1f%% of the ',type,' support matches ',atitle0,'\n'],length(iy)/length(y0)*100)
- icheck = mcf_analysis(x0, iy, iyc, options_mcf, options_, bayestopt_, estim_params_);
+ icheck = mcf_analysis(x0, iy, iyc, options_mcf, M_, options_, bayestopt_, estim_params_);
lpmat=x0(iy,:);
if nshocks
diff --git a/matlab/gsa/redform_screen.m b/matlab/gsa/redform_screen.m
index 29a7f3c525e7645722971c3a76b4b7056d5a399c..40dd532879d04acbf7c856b5738a0f6655e81f22 100644
--- a/matlab/gsa/redform_screen.m
+++ b/matlab/gsa/redform_screen.m
@@ -1,10 +1,10 @@
-function redform_screen(dirname, options_gsa_, param_vals, M_, dr, options_, bayestopt_)
-% redform_screen(dirname, options_gsa_, param_vals, M_, dr, options_, bayestopt_)
+function redform_screen(dirname, options_gsa_, estim_params_, M_, dr, options_, bayestopt_)
+% redform_screen(dirname, options_gsa_, estim_params_, M_, dr, options_, bayestopt_)
% Conduct reduced form screening
% Inputs:
% - dirname [string] name of the output directory
% - options_gsa_ [structure] GSA options_
-% - param_vals [double] parameter vector
+% - estim_params [structure] describing the estimated parameters
% - M_ [structure] describing the model
% - dr [structure] decision rules
% - options_ [structure] describing the options
@@ -37,7 +37,12 @@ anamlagendo = options_gsa_.namlagendo;
anamexo = options_gsa_.namexo;
nliv = options_gsa_.morris_nliv;
-pnames = M_.param_names(param_vals(:,1));
+pnames = M_.param_names(estim_params_.param_vals(:,1));
+if options_.TeX
+ for par_iter=1:size(estim_params_.param_vals(:,1),1)
+ [~,tex_names{par_iter,1}]=get_the_name(estim_params_.param_vals(par_iter,1),options_.TeX, M_, estim_params_, options_);
+ end
+end
if nargin==0
dirname='';
end
@@ -53,12 +58,12 @@ nsok = length(find(M_.lead_lag_incidence(M_.maximum_lag,:)));
js=0;
for j=1:size(anamendo,1)
- namendo = deblank(anamendo(j,:));
+ namendo = anamendo{j,:};
iendo = strmatch(namendo, M_.endo_names(dr.order_var), 'exact');
iplo=0;
ifig=0;
for jx=1:size(anamexo,1)
- namexo = deblank(anamexo(jx,:));
+ namexo = anamexo{jx};
iexo = strmatch(namexo, M_.exo_names, 'exact');
if ~isempty(iexo)
y0=teff(T(iendo,iexo+nspred,:), kn, istable);
@@ -79,7 +84,11 @@ for j=1:size(anamendo,1)
set(gca,'xticklabel',' ','fontsize',10)
set(gca,'xlim',[0.5 10.5])
for ip=1:min(np,10)
- text(ip,-0.02,pnames(iso(ip)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
+ if options_.TeX
+ text(ip,-0.02,tex_names(iso(ip)),'rotation',90,'HorizontalAlignment','right','interpreter','latex')
+ else
+ text(ip,-0.02,pnames(iso(ip)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
+ end
end
title([namendo,' vs. ',namexo],'interpreter','none')
if iplo==9
@@ -98,7 +107,7 @@ for j=1:size(anamendo,1)
iplo=0;
ifig=0;
for je=1:size(anamlagendo,1)
- namlagendo=deblank(anamlagendo(je,:));
+ namlagendo=anamlagendo{je};
ilagendo=strmatch(namlagendo, M_.endo_names(dr.order_var(M_.nstatic+1:M_.nstatic+nsok)), 'exact');
if ~isempty(ilagendo)
@@ -120,7 +129,11 @@ for j=1:size(anamendo,1)
set(gca,'xticklabel',' ','fontsize',10)
set(gca,'xlim',[0.5 10.5])
for ip=1:min(np,10)
- text(ip,-0.02,deblank(pnames(iso(ip),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
+ if options_.TeX
+ text(ip,-0.02,tex_names(iso(ip)),'rotation',90,'HorizontalAlignment','right','interpreter','latex')
+ else
+ text(ip,-0.02,pnames{iso(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
+ end
end
title([namendo,' vs. ',namlagendo,'(-1)'],'interpreter','none')
@@ -144,7 +157,11 @@ set(gca,'xlim',[0.5 np+0.5])
set(gca,'ylim',[0 1])
set(gca,'position',[0.13 0.2 0.775 0.7])
for ip=1:np
- text(ip,-0.02,deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
+ if options_.TeX
+ text(ip,-0.02,tex_names(ip),'rotation',90,'HorizontalAlignment','right','interpreter','latex')
+ else
+ text(ip,-0.02,pnames{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
+ end
end
xlabel(' ')
ylabel('Elementary Effects')
diff --git a/matlab/gsa/scatter_mcf.m b/matlab/gsa/scatter_mcf.m
index fbb40dd24da77c4ed8ca14a7672f99fc3f805630..8c6c1908f44466d01552fe87cf35df0c4300459f 100644
--- a/matlab/gsa/scatter_mcf.m
+++ b/matlab/gsa/scatter_mcf.m
@@ -38,12 +38,6 @@ function scatter_mcf(X,Y,vnames,plotsymbol, fnam, dirname, figtitle, xparam1, o
Z=[X;Y];
-[n,p] = size(X);
-% X = X - ones(n,1)*min(Z);
-% X = X ./ (ones(n,1)*max(Z));
-[n,p] = size(Y);
-% Y = Y - ones(n,1)*min(Z);
-% Y = Y ./ (ones(n,1)*max(Z));
[n,p] = size(Z);
clear Z;
@@ -53,8 +47,10 @@ if nargin >=3
end
if nargin<4 || isempty(plotsymbol)
- if n*p<100, plotsymbol = 'o';
- else plotsymbol = '.';
+ if n*p<100
+ plotsymbol = 'o';
+ else
+ plotsymbol = '.';
end
end
@@ -83,7 +79,7 @@ end
figtitle_tex=strrep(figtitle,'_','\_');
-fig_nam_=[fnam];
+fig_nam_=fnam;
if ~nograph
hh_fig=dyn_figure(options_.nodisplay,'name',figtitle);
end
@@ -101,7 +97,6 @@ for i = 1:p
h = axes('position',[fL(i),fL(p+1-j),ffl,ffl]);
if i==j
h1=cumplot(X(:,j));
- % set(h1,'color',[0 0 1], 'linestyle','--','LineWidth',1.5)
set(h1,'color',[0 0 1],'LineWidth',1.5)
hold on,
h2=cumplot(Y(:,j));
@@ -126,10 +121,10 @@ for i = 1:p
plot(X(:,i),X(:,j),[plotsymbol,'b'])
end
if ~isempty(xparam1)
- hold on, plot(xparam1(i),xparam1(j),'s','MarkerFaceColor',[0 0.75 0],'MarkerEdgeColor',[0 0.75 0])
+ hold on
+ plot(xparam1(i),xparam1(j),'s','MarkerFaceColor',[0 0.75 0],'MarkerEdgeColor',[0 0.75 0])
end
hold off;
- % axis([-0.1 1.1 -0.1 1.1])
if i<p
set(gca,'YTickLabel',[],'YTick',[]);
else
@@ -144,16 +139,26 @@ for i = 1:p
end
if i==1
if nflag == 1
- ylabel(vnames(j,:),'Rotation',45, ...
- 'HorizontalAlignment','right','VerticalAlignment','middle');
+ if options_.TeX
+ ylabel(vnames(j,:),'Rotation',45, ...
+ 'HorizontalAlignment','right','VerticalAlignment','middle','Interpreter','latex');
+ else
+ ylabel(vnames(j,:),'Rotation',45, ...
+ 'HorizontalAlignment','right','VerticalAlignment','middle','Interpreter','none');
+ end
else
ylabel([num2str(j),' '],'Rotation',90)
end
end
if j==1
if nflag == 1
- title(vnames(i,:),'Rotation',45, ...
- 'HorizontalAlignment','left','VerticalAlignment','bottom')
+ if options_.TeX
+ title(vnames(i,:),'Rotation',45, ...
+ 'HorizontalAlignment','left','VerticalAlignment','bottom','Interpreter','latex')
+ else
+ title(vnames(i,:),'Rotation',45, ...
+ 'HorizontalAlignment','left','VerticalAlignment','bottom','Interpreter','none')
+ end
else
title(num2str(i))
end
@@ -162,8 +167,13 @@ for i = 1:p
end
end
if ~isoctave
- annotation('textbox', [0.1,0,0.35,0.05],'String', beha_name,'Color','Blue','horizontalalignment','center','interpreter','none');
- annotation('textbox', [0.55,0,0.35,0.05],'String', non_beha_name,'Color','Red','horizontalalignment','center','interpreter','none');
+ if options_.TeX
+ annotation('textbox', [0.1,0,0.35,0.05],'String', beha_name,'Color','Blue','horizontalalignment','center','interpreter','latex');
+ annotation('textbox', [0.55,0,0.35,0.05],'String', non_beha_name,'Color','Red','horizontalalignment','center','interpreter','latex');
+ else
+ annotation('textbox', [0.1,0,0.35,0.05],'String', beha_name,'Color','Blue','horizontalalignment','center','interpreter','none');
+ annotation('textbox', [0.55,0,0.35,0.05],'String', non_beha_name,'Color','Red','horizontalalignment','center','interpreter','none');
+ end
end
if ~nograph
diff --git a/matlab/gsa/set_shocks_param.m b/matlab/gsa/set_shocks_param.m
index 66141038c480d860b527ba163998e51cb7925479..02d85630b4814ea1325b2b55a222b5661b41eb4e 100644
--- a/matlab/gsa/set_shocks_param.m
+++ b/matlab/gsa/set_shocks_param.m
@@ -1,8 +1,16 @@
-function set_shocks_param(xparam1)
-% function set_shocks_param(xparam1)
+function M_=set_shocks_param(M_,estim_params_,xparam1)
+% function M_=set_shocks_param(M_,estim_params_,xparam1)
% Set the structural and measurement error variances and covariances
+% Inputs
+% - M_ [structure] Matlab's structure describing the model
+% - estim_params_ [structure] characterizing parameters to be estimated
+% - xparam1 [double] parameter vector
+% Outputs:
+% - M_ [structure] Matlab's structure describing the model
+%
+% Notes: closely follows set_all_parameters.m
-% Copyright © 2012-2017 Dynare Team
+% Copyright © 2012-2023 Dynare Team
%
% This file is part of Dynare.
%
@@ -19,8 +27,6 @@ function set_shocks_param(xparam1)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-global estim_params_ M_
-
nvx = estim_params_.nvx;
ncx = estim_params_.ncx;
nvn = estim_params_.nvn;
diff --git a/matlab/gsa/stab_map_.m b/matlab/gsa/stab_map_.m
index 3dc6d1082918ac59539c933f84d1c334bc892230..53521dac19116f2d6fd1cda28505872c3857c8f5 100644
--- a/matlab/gsa/stab_map_.m
+++ b/matlab/gsa/stab_map_.m
@@ -1,9 +1,22 @@
-function x0 = stab_map_(OutputDirectoryName,opt_gsa)
-% x0 = stab_map_(OutputDirectoryName,opt_gsa)
+function x0 = stab_map_(OutputDirectoryName,opt_gsa,M_,oo_,options_,bayestopt_,estim_params_)
+% x0 = stab_map_(OutputDirectoryName,opt_gsa,M_,oo_,options_,bayestopt_,estim_params_)
% Mapping of stability regions in the prior ranges applying
% Monte Carlo filtering techniques.
%
-% INPUTS (from opt_gsa structure)
+% Inputs
+% - OutputDirectoryName [string] name of the output directory
+% - opt_gsa [structure] GSA options structure
+% - M_ [structure] Matlab's structure describing the model
+% - oo_ [structure] Matlab's structure describing the results
+% - options_ [structure] Matlab's structure describing the current options
+% - bayestopt_ [structure] describing the priors
+% - estim_params_ [structure] characterizing parameters to be estimated
+%
+% Outputs:
+% - x0 one parameter vector for which the model is stable.
+%
+%
+% Inputs from opt_gsa structure
% Nsam = MC sample size
% fload = 0 to run new MC; 1 to load prevoiusly generated analysis
% alpha2 = significance level for bivariate sensitivity analysis
@@ -14,8 +27,6 @@ function x0 = stab_map_(OutputDirectoryName,opt_gsa)
% _prior.mat file
% = 0: sample from posterior ranges: sample saved in
% _mc.mat file
-% OUTPUT:
-% x0: one parameter vector for which the model is stable.
%
% GRAPHS
% 1) Pdf's of marginal distributions under the stability (dotted
@@ -50,11 +61,6 @@ function x0 = stab_map_(OutputDirectoryName,opt_gsa)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-%global bayestopt_ estim_params_ dr_ options_ ys_ fname_
-global bayestopt_ estim_params_ options_ oo_ M_
-
-% opt_gsa=options_.opt_gsa;
-
Nsam = opt_gsa.Nsam;
fload = opt_gsa.load_stab;
alpha2 = opt_gsa.alpha2_stab;
@@ -82,6 +88,7 @@ nshock = nshock + estim_params_.ncn;
lpmat0=zeros(Nsam,0);
xparam1=[];
+%% prepare prior bounds
[~,~,~,lb,ub,~] = set_prior(estim_params_,M_,options_); %Prepare bounds
if ~isempty(bayestopt_) && any(bayestopt_.pshape > 0)
% Set prior bounds
@@ -107,12 +114,13 @@ options_mcf.pvalue_ks = pvalue_ks;
options_mcf.pvalue_corr = pvalue_corr;
options_mcf.alpha2 = alpha2;
+%% get LaTeX names
name=cell(np,1);
name_tex=cell(np,1);
for jj=1:np
if options_.TeX
[param_name_temp, param_name_tex_temp]= get_the_name(nshock+jj,options_.TeX,M_,estim_params_,options_);
- name_tex{jj,1} = strrep(param_name_tex_temp,'$','');
+ name_tex{jj,1} = param_name_tex_temp;
name{jj,1} = param_name_temp;
else
param_name_temp = get_the_name(nshock+jj,options_.TeX,M_,estim_params_,options_);
@@ -128,19 +136,18 @@ options_mcf.fname_ = fname_;
options_mcf.OutputDirectoryName = OutputDirectoryName;
options_mcf.xparam1 = [];
-opt=options_;
options_.periods=0;
options_.nomoments=1;
options_.irf=0;
options_.noprint=1;
-if fload==0
+if fload==0 %run new MC
if isfield(dr_,'ghx')
egg=zeros(length(dr_.eigval),Nsam);
end
yys=zeros(length(dr_.ys),Nsam);
if opt_gsa.morris == 1
- [lpmat, OutFact] = Sampling_Function_2(nliv, np+nshock, ntra, ones(np+nshock, 1), zeros(np+nshock,1), []);
+ [lpmat] = Sampling_Function_2(nliv, np+nshock, ntra, ones(np+nshock, 1), zeros(np+nshock,1), []);
lpmat = lpmat.*(nliv-1)/nliv+1/nliv/2;
Nsam=size(lpmat,1);
lpmat0 = lpmat(:,1:nshock);
@@ -158,10 +165,9 @@ if fload==0
for j=1:np
lpmat(:,j) = randperm(Nsam)'./(Nsam+1); %latin hypercube
end
-
end
end
- dummy=prior_draw_gsa(1);
+ prior_draw_gsa(M_,bayestopt_,options_,estim_params_,1); %initialize
if pprior
for j=1:nshock
if opt_gsa.morris~=1
@@ -178,16 +184,16 @@ if fload==0
lpmat(:,j)=lpmat(:,j).*(upper_bound-lower_bound)+lower_bound;
end
else
- xx=prior_draw_gsa(0,[lpmat0 lpmat]);
+ xx=prior_draw_gsa(M_,bayestopt_,options_,estim_params_,0,[lpmat0 lpmat]);
lpmat0=xx(:,1:nshock);
lpmat=xx(:,nshock+1:end);
clear xx;
end
- else
+ else %posterior analysis
if neighborhood_width>0 && isempty(options_.mode_file)
xparam1 = get_all_parameters(estim_params_,M_);
else
- eval(['load ' options_.mode_file '.mat;']);
+ load([options_.mode_file '.mat'],'hh','xparam1');
end
if neighborhood_width>0
for j=1:nshock
@@ -216,8 +222,8 @@ if fload==0
for j=1:Nsam*2
lnprior(j) = any(lp(j,:)'<=bounds.lb | lp(j,:)'>=bounds.ub);
end
- ireal=[1:2*Nsam];
- ireal=ireal(find(lnprior==0));
+ ireal=1:2*Nsam;
+ ireal=ireal(lnprior==0);
lp=lp(ireal,:);
Nsam=min(Nsam, length(ireal));
lpmat0=lp(1:Nsam,1:nshock);
@@ -227,9 +233,9 @@ if fload==0
end
%
h = dyn_waitbar(0,'Please wait...');
- istable=[1:Nsam];
+ istable=1:Nsam;
jstab=0;
- iunstable=[1:Nsam];
+ iunstable=1:Nsam;
iindeterm=zeros(1,Nsam);
iwrong=zeros(1,Nsam);
inorestriction=zeros(1,Nsam);
@@ -237,12 +243,11 @@ if fload==0
infox=zeros(Nsam,1);
for j=1:Nsam
M_ = set_all_parameters([lpmat0(j,:) lpmat(j,:)]',estim_params_,M_);
- %try stoch_simul([]);
try
- if ~ isempty(options_.endogenous_prior_restrictions.moment)
- [Tt,Rr,SteadyState,info,oo_.dr,M_.params] = dynare_resolve(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
+ if ~isempty(options_.endogenous_prior_restrictions.moment)
+ [Tt,Rr,~,info,oo_.dr,M_.params] = dynare_resolve(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
else
- [Tt,Rr,SteadyState,info,oo_.dr,M_.params] = dynare_resolve(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state,'restrict');
+ [Tt,Rr,~,info,oo_.dr,M_.params] = dynare_resolve(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state,'restrict');
end
infox(j,1)=info(1);
if infox(j,1)==0 && ~exist('T','var')
@@ -250,8 +255,7 @@ if fload==0
if prepSA
try
T=zeros(size(dr_.ghx,1),size(dr_.ghx,2)+size(dr_.ghu,2),Nsam);
- catch
- ME = lasterror();
+ catch ME
if strcmp('MATLAB:nomem',ME.identifier)
prepSA=0;
disp('The model is too large for storing state space matrices ...')
@@ -324,7 +328,6 @@ if fload==0
end
ys_=real(dr_.ys);
yys(:,j) = ys_;
- ys_=yys(:,1);
dyn_waitbar(j/Nsam,h,['MC iteration ',int2str(j),'/',int2str(Nsam)])
end
dyn_waitbar_close(h);
@@ -333,13 +336,13 @@ if fload==0
else
T=[];
end
- istable=istable(find(istable)); % stable params ignoring restrictions
- irestriction=irestriction(find(irestriction)); % stable params & restrictions OK
- inorestriction=inorestriction(find(inorestriction)); % stable params violating restrictions
- iunstable=iunstable(find(iunstable)); % violation of BK & restrictions & solution could not be found (whatever goes wrong)
- iindeterm=iindeterm(find(iindeterm)); % indeterminacy
- iwrong=iwrong(find(iwrong)); % dynare could not find solution
- ixun=iunstable(find(~ismember(iunstable,[iindeterm,iwrong,inorestriction]))); % explosive roots
+ istable=istable(istable~=0); % stable params ignoring restrictions
+ irestriction=irestriction(irestriction~=0); % stable params & restrictions OK
+ inorestriction=inorestriction(inorestriction~=0); % stable params violating restrictions
+ iunstable=iunstable(iunstable~=0); % violation of BK & restrictions & solution could not be found (whatever goes wrong)
+ iindeterm=iindeterm(iindeterm~=0); % indeterminacy
+ iwrong=iwrong(iwrong~=0); % dynare could not find solution
+ ixun=iunstable(~ismember(iunstable,[iindeterm,iwrong,inorestriction])); % explosive roots
bkpprior.pshape=bayestopt_.pshape;
bkpprior.p1=bayestopt_.p1;
@@ -356,8 +359,7 @@ if fload==0
'bkpprior','lpmat','lpmat0','irestriction','iunstable','istable','iindeterm','iwrong','ixun', ...
'egg','yys','T','nspred','nboth','nfwrd','infox')
end
-
- else
+ else %~pprior
if ~prepSA
save([OutputDirectoryName filesep fname_ '_mc.mat'], ...
'lpmat','lpmat0','irestriction','iunstable','istable','iindeterm','iwrong','ixun', ...
@@ -368,19 +370,18 @@ if fload==0
'egg','yys','T','nspred','nboth','nfwrd','infox')
end
end
-else
+else %load old run
if pprior
filetoload=[OutputDirectoryName filesep fname_ '_prior.mat'];
else
filetoload=[OutputDirectoryName filesep fname_ '_mc.mat'];
end
- load(filetoload,'lpmat','lpmat0','irestriction','iunstable','istable','iindeterm','iwrong','ixun','egg','yys','nspred','nboth','nfwrd','infox')
+ load(filetoload,'lpmat','lpmat0','irestriction','iunstable','istable','iindeterm','iwrong','ixun','infox')
Nsam = size(lpmat,1);
if pprior==0 && ~isempty(options_.mode_file)
- eval(['load ' options_.mode_file '.mat;']);
+ load([options_.mode_file '.mat'],'xparam1');
end
-
if prepSA && isempty(strmatch('T',who('-file', filetoload),'exact'))
h = dyn_waitbar(0,'Please wait...');
options_.periods=0;
@@ -392,30 +393,22 @@ else
yys=NaN(length(ys_),ntrans);
for j=1:ntrans
M_.params(estim_params_.param_vals(:,1)) = lpmat(istable(j),:)';
- %stoch_simul([]);
- [Tt,Rr,SteadyState,info,oo_.dr,M_.params] = dynare_resolve(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state,'restrict');
+ [~,~,~,~,oo_.dr,M_.params] = dynare_resolve(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state,'restrict');
if ~exist('T','var')
T=zeros(size(dr_.ghx,1),size(dr_.ghx,2)+size(dr_.ghu,2),ntrans);
end
dr_ = oo_.dr;
T(:,:,j) = [dr_.ghx dr_.ghu];
- if ~exist('nspred','var')
- nspred = dr_.nspred; %size(dr_.ghx,2);
- nboth = dr_.nboth;
- nfwrd = dr_.nfwrd;
- end
ys_=real(dr_.ys);
yys(:,j) = ys_;
- ys_=yys(:,1);
dyn_waitbar(j/ntrans,h,['MC iteration ',int2str(j),'/',int2str(ntrans)])
end
dyn_waitbar_close(h);
save(filetoload,'T','-append')
- elseif prepSA
- load(filetoload,'T')
end
end
+%% display and save output
if pprior
aunstname='prior_unstable'; aunsttitle='Prior StabMap: explosiveness of solution';
aindname='prior_indeterm'; aindtitle='Prior StabMap: Indeterminacy';
@@ -435,17 +428,18 @@ delete([OutputDirectoryName,filesep,fname_,'_',aindname,'.*']);
delete([OutputDirectoryName,filesep,fname_,'_',aunstname,'.*']);
delete([OutputDirectoryName,filesep,fname_,'_',awrongname,'.*']);
-if length(iunstable)>0 || length(iwrong)>0
- fprintf(['%4.1f%% of the prior support gives unique saddle-path solution.\n'],length(istable)/Nsam*100)
- fprintf(['%4.1f%% of the prior support gives explosive dynamics.\n'],(length(ixun) )/Nsam*100)
+fprintf('\nSensitivity Analysis: Stability mapping:\n')
+if ~isempty(iunstable) || ~isempty(iwrong)
+ fprintf('%4.1f%% of the prior support gives unique saddle-path solution.\n',length(istable)/Nsam*100)
+ fprintf('%4.1f%% of the prior support gives explosive dynamics.\n',(length(ixun) )/Nsam*100)
if ~isempty(iindeterm)
- fprintf(['%4.1f%% of the prior support gives indeterminacy.\n'],length(iindeterm)/Nsam*100)
+ fprintf('%4.1f%% of the prior support gives indeterminacy.\n',length(iindeterm)/Nsam*100)
end
- inorestriction = istable(find(~ismember(istable,irestriction))); % violation of prior restrictions
+ inorestriction = istable(~ismember(istable,irestriction)); % violation of prior restrictions
if ~isempty(iwrong) || ~isempty(inorestriction)
skipline()
if any(infox==49)
- fprintf(['%4.1f%% of the prior support violates prior restrictions.\n'],(length(inorestriction) )/Nsam*100)
+ fprintf('%4.1f%% of the prior support violates prior restrictions.\n',(length(inorestriction) )/Nsam*100)
end
if ~isempty(iwrong)
skipline()
@@ -486,50 +480,50 @@ if length(iunstable)>0 || length(iwrong)>0
end
skipline()
if length(iunstable)<Nsam || length(istable)>1
- itot = [1:Nsam];
- isolve = itot(find(~ismember(itot,iwrong))); % dynare could find a solution
+ itot = 1:Nsam;
+ isolve = itot(~ismember(itot,iwrong)); % dynare could find a solution
% Blanchard Kahn
if neighborhood_width
options_mcf.xparam1 = xparam1(nshock+1:end);
end
- itmp = itot(find(~ismember(itot,istable)));
+ itmp = itot(~ismember(itot,istable));
options_mcf.amcf_name = asname;
options_mcf.amcf_title = atitle;
options_mcf.beha_title = 'unique Stable Saddle-Path';
options_mcf.nobeha_title = 'NO unique Stable Saddle-Path';
options_mcf.title = 'unique solution';
- mcf_analysis(lpmat, istable, itmp, options_mcf, options_, bayestopt_, estim_params_);
+ mcf_analysis(lpmat, istable, itmp, options_mcf, M_, options_, bayestopt_, estim_params_);
if ~isempty(iindeterm)
- itmp = isolve(find(~ismember(isolve,iindeterm)));
+ itmp = isolve(~ismember(isolve,iindeterm));
options_mcf.amcf_name = aindname;
options_mcf.amcf_title = aindtitle;
options_mcf.beha_title = 'NO indeterminacy';
options_mcf.nobeha_title = 'indeterminacy';
options_mcf.title = 'indeterminacy';
- mcf_analysis(lpmat, itmp, iindeterm, options_mcf, options_, bayestopt_, estim_params_);
+ mcf_analysis(lpmat, itmp, iindeterm, options_mcf, M_, options_, bayestopt_, estim_params_);
end
if ~isempty(ixun)
- itmp = isolve(find(~ismember(isolve,ixun)));
+ itmp = isolve(~ismember(isolve,ixun));
options_mcf.amcf_name = aunstname;
options_mcf.amcf_title = aunsttitle;
options_mcf.beha_title = 'NO explosive solution';
options_mcf.nobeha_title = 'explosive solution';
options_mcf.title = 'instability';
- mcf_analysis(lpmat, itmp, ixun, options_mcf, options_, bayestopt_, estim_params_);
+ mcf_analysis(lpmat, itmp, ixun, options_mcf, M_, options_, bayestopt_, estim_params_);
end
- inorestriction = istable(find(~ismember(istable,irestriction))); % violation of prior restrictions
- iwrong = iwrong(find(~ismember(iwrong,inorestriction))); % what went wrong beyond prior restrictions
+ inorestriction = istable(~ismember(istable,irestriction)); % violation of prior restrictions
+ iwrong = iwrong(~ismember(iwrong,inorestriction)); % what went wrong beyond prior restrictions
if ~isempty(iwrong)
- itmp = itot(find(~ismember(itot,iwrong)));
+ itmp = itot(~ismember(itot,iwrong));
options_mcf.amcf_name = awrongname;
options_mcf.amcf_title = awrongtitle;
options_mcf.beha_title = 'NO inability to find a solution';
options_mcf.nobeha_title = 'inability to find a solution';
options_mcf.title = 'inability to find a solution';
- mcf_analysis(lpmat, itmp, iwrong, options_mcf, options_, bayestopt_, estim_params_);
+ mcf_analysis(lpmat, itmp, iwrong, options_mcf, M_, options_, bayestopt_, estim_params_);
end
if ~isempty(irestriction)
@@ -542,7 +536,7 @@ if length(iunstable)>0 || length(iwrong)>0
for jj=1:np
if options_.TeX
[param_name_temp, param_name_tex_temp]= get_the_name(jj,options_.TeX,M_,estim_params_,options_);
- name_tex{jj,1} = strrep(param_name_tex_temp,'$','');
+ name_tex{jj,1} = param_name_tex_temp;
name{jj,1} = param_name_temp;
else
param_name_temp = get_the_name(jj,options_.TeX,M_,estim_params_,options_);
@@ -558,7 +552,7 @@ if length(iunstable)>0 || length(iwrong)>0
options_mcf.beha_title = 'prior IRF/moment calibration';
options_mcf.nobeha_title = 'NO prior IRF/moment calibration';
options_mcf.title = 'prior restrictions';
- mcf_analysis([lpmat0 lpmat], irestriction, inorestriction, options_mcf, options_, bayestopt_, estim_params_);
+ mcf_analysis([lpmat0 lpmat], irestriction, inorestriction, options_mcf, M_, options_, bayestopt_, estim_params_);
iok = irestriction(1);
x0 = [lpmat0(iok,:)'; lpmat(iok,:)'];
else
@@ -568,7 +562,7 @@ if length(iunstable)>0 || length(iwrong)>0
end
M_ = set_all_parameters(x0,estim_params_,M_);
- [oo_.dr,info,M_.params] = resol(0,M_,options_,oo_.dr ,oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
+ [oo_.dr,~,M_.params] = resol(0,M_,options_,oo_.dr ,oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
else
disp('All parameter values in the specified ranges are not acceptable!')
x0=[];
@@ -578,15 +572,8 @@ else
disp('and match prior IRF/moment restriction(s) if any!')
x0=0.5.*(bounds.ub(1:nshock)-bounds.lb(1:nshock))+bounds.lb(1:nshock);
x0 = [x0; lpmat(istable(1),:)'];
-
end
+skipline(1);
xparam1=x0;
-save([OutputDirectoryName filesep 'prior_ok.mat'],'xparam1');
-
-options_.periods=opt.periods;
-if isfield(opt,'nomoments')
- options_.nomoments=opt.nomoments;
-end
-options_.irf=opt.irf;
-options_.noprint=opt.noprint;
+save([OutputDirectoryName filesep 'prior_ok.mat'],'xparam1');
\ No newline at end of file
diff --git a/matlab/gsa/stab_map_2.m b/matlab/gsa/stab_map_2.m
index 50a33603f6725787cc03b9449699967f55245ca2..f7d0ca8bd64ac052eeff70acb80da2863ec73456 100644
--- a/matlab/gsa/stab_map_2.m
+++ b/matlab/gsa/stab_map_2.m
@@ -1,6 +1,20 @@
-function indcorr = stab_map_2(x,alpha2, pvalue_crit, fnam, dirname,xparam1,figtitle)
-% function stab_map_2(x, alpha2, pvalue, fnam, dirname,xparam1)
+function indcorr = stab_map_2(x,alpha2, pvalue_crit, M_,options_,bayestopt_,estim_params_, fnam, dirname,xparam1,figtitle)
+% indcorr = stab_map_2(x,alpha2, pvalue_crit, M_,options_,bayestopt_,estim_params_, fnam, dirname,xparam1,figtitle)
+% Inputs:
+% - x
+% - alpha2
+% - pvalue_crit
+% - M_ [structure] Matlab's structure describing the model
+% - options_ [structure] Matlab's structure describing the current options
+% - bayestopt_ [structure] describing the priors
+% - estim_params_ [structure] characterizing parameters to be estimated
+% - fnam [string] file name
+% - dirnam [string] directory name
+% - xparam1 [double] parameter vector
+% - figtitle [string] figure title
%
+% Output:
+% - indcorr
% Written by Marco Ratto
% Joint Research Centre, The European Commission,
% marco.ratto@ec.europa.eu
@@ -22,29 +36,23 @@ function indcorr = stab_map_2(x,alpha2, pvalue_crit, fnam, dirname,xparam1,figti
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-%global bayestopt_ estim_params_ dr_ options_ ys_ fname_
-global bayestopt_ estim_params_ options_ oo_ M_
-
npar=size(x,2);
-nsam=size(x,1);
ishock= npar>estim_params_.np;
nograph = options_.nograph;
-if nargin<4
+if nargin<8
fnam='';
end
-if nargin<5
+if nargin<9
dirname='';
nograph=1;
end
-if nargin<6
+if nargin<10
xparam1=[];
end
-if nargin<7
+if nargin<11
figtitle=fnam;
end
-ys_ = oo_.dr.ys;
-dr_ = oo_.dr;
fname_ = M_.fname;
nshock = estim_params_.nvx;
nshock = nshock + estim_params_.nvn;
@@ -74,7 +82,7 @@ indcorr = [];
entry_iter=1;
for j=1:npar
i2=find(abs(c00(:,j))>alpha2);
- if length(i2)>0
+ if ~isempty(i2)
for jx=1:length(i2)
if pvalue(j,i2(jx))<pvalue_crit
indcorr = [indcorr; [j i2(jx)]];
@@ -82,9 +90,7 @@ for j=1:npar
if ishock
if options_.TeX
[param_name_temp1, param_name_tex_temp1]= get_the_name(j,options_.TeX,M_,estim_params_,options_);
- param_name_tex_temp1 = strrep(param_name_tex_temp1,'$','');
[param_name_temp2, param_name_tex_temp2]= get_the_name(i2(jx),options_.TeX,M_,estim_params_,options_);
- param_name_tex_temp2 = strrep(param_name_tex_temp2,'$','');
tmp_name=(['[',param_name_temp1,',',param_name_temp2,']']);
tmp_name_tex=(['[',param_name_tex_temp1,',',param_name_tex_temp2,']']);
name{entry_iter,1}=tmp_name;
@@ -98,9 +104,7 @@ for j=1:npar
else
if options_.TeX
[param_name_temp1, param_name_tex_temp1]= get_the_name(j+nshock,options_.TeX,M_,estim_params_,options_);
- param_name_tex_temp1 = strrep(param_name_tex_temp1,'$','');
[param_name_temp2, param_name_tex_temp2]= get_the_name(i2(jx)+nshock,options_.TeX,M_,estim_params_,options_);
- param_name_tex_temp2 = strrep(param_name_tex_temp2,'$','');
tmp_name=(['[',param_name_temp1,',',param_name_temp2,']']);
tmp_name_tex=(['[',param_name_tex_temp1,',',param_name_tex_temp2,']']);
name{entry_iter,1}=tmp_name;
@@ -121,17 +125,10 @@ for j=1:npar
hh_fig=dyn_figure(options_.nodisplay,'name',[figtitle,' sample bivariate projection ', num2str(ifig)]);
end
subplot(3,4,j2-(ifig-1)*12)
- % bar(c0(i2,j)),
- % set(gca,'xticklabel',bayestopt_.name(i2)),
- % set(gca,'xtick',[1:length(i2)])
- %plot(stock_par(ixx(nfilt+1:end,i),j),stock_par(ixx(nfilt+1:end,i),i2(jx)),'.k')
- %hold on,
plot(x(:,j),x(:,i2(jx)),'.')
if ~isempty(xparam1)
hold on, plot(xparam1(j),xparam1(i2(jx)),'ro')
end
- % xlabel(deblank(estim_params_.param_names(j,:)),'interpreter','none'),
- % ylabel(deblank(estim_params_.param_names(i2(jx),:)),'interpreter','none'),
if ishock
xlabel(bayestopt_.name{j},'interpreter','none'),
ylabel(bayestopt_.name{i2(jx)},'interpreter','none'),
@@ -191,5 +188,4 @@ else
if options_.TeX
dyn_latex_table(M_, options_, title_string_tex, fig_nam_tex_table, headers, name_tex, data_mat, 0, 7, 3);
end
-end
-%close all
+end
\ No newline at end of file
diff --git a/matlab/plot_identification.m b/matlab/plot_identification.m
index f0d50e4c592cc8e9dded052b04bf325ef5a5ee35..fcfbc60748f3b28678009d6e14ed0bce4ea551d9 100644
--- a/matlab/plot_identification.m
+++ b/matlab/plot_identification.m
@@ -1,7 +1,8 @@
-function plot_identification(params, idemoments, idehess, idemodel, idelre, advanced, tittxt, name, IdentifDirectoryName, fname, options_, estim_params_, tit_TeX, name_tex)
-% function plot_identification(params,idemoments,idehess,idemodel, idelre, advanced, tittxt, name, IdentifDirectoryName)
+function plot_identification(M_, params, idemoments, idehess, idemodel, idelre, advanced, tittxt, name, IdentifDirectoryName, fname, options_, estim_params_, bayestopt_, tit_TeX, name_tex)
+% function plot_identification(params,idemoments,idehess,idemodel, idelre, advanced, tittxt, name, IdentifDirectoryName, fname, options_, estim_params_, bayestopt_, tit_TeX, name_tex)
%
% INPUTS
+% o M_ [structure] model
% o params [array] parameter values for identification checks
% o idemoments [structure] identification results for the moments
% o idehess [structure] identification results for the Hessian
@@ -36,11 +37,11 @@ function plot_identification(params, idemoments, idehess, idemodel, idelre, adva
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-if nargin <12 || isempty(tit_TeX)
+if nargin <14 || isempty(tit_TeX)
tit_TeX=tittxt;
end
-if nargin <13
+if nargin <15
name_tex=name;
end
@@ -436,17 +437,17 @@ else
options_mcf.title = 'MC Highest Condition Number LRE Model';
ncut=floor(SampleSize/10*9);
[~,is]=sort(idelre.cond);
- mcf_analysis(params, is(1:ncut), is(ncut+1:end), options_mcf, options_, bayestopt_, estim_params_);
+ mcf_analysis(params, is(1:ncut), is(ncut+1:end), options_mcf, M_, options_, bayestopt_, estim_params_);
options_mcf.amcf_name = 'MC_HighestCondNumberModel';
options_mcf.amcf_title = 'MC Highest Condition Number Model Solution';
options_mcf.title = 'MC Highest Condition Number Model Solution';
[~,is]=sort(idemodel.cond);
- mcf_analysis(params, is(1:ncut), is(ncut+1:end), options_mcf, options_, bayestopt_, estim_params_);
+ mcf_analysis(params, is(1:ncut), is(ncut+1:end), options_mcf, M_, options_, bayestopt_, estim_params_);
options_mcf.amcf_name = 'MC_HighestCondNumberMoments';
options_mcf.amcf_title = 'MC Highest Condition Number Model Moments';
options_mcf.title = 'MC Highest Condition Number Model Moments';
[~,is]=sort(idemoments.cond);
- mcf_analysis(params, is(1:ncut), is(ncut+1:end), options_mcf, options_, bayestopt_, estim_params_);
+ mcf_analysis(params, is(1:ncut), is(ncut+1:end), options_mcf, M_, options_, bayestopt_, estim_params_);
% [proba, dproba] = stab_map_1(idemoments.Mco', is(1:ncut), is(ncut+1:end), 'HighestCondNumberMoments_vs_Mco', 1, [], IdentifDirectoryName);
% for j=1:nparam,
% % ibeh=find(idemoments.Mco(j,:)<0.9);
diff --git a/preprocessor b/preprocessor
index 6d9fc367d034be8db19857995f413fd8a740ebba..378d00fc3a50a11c3dce615ea337ad55989c938e 160000
--- a/preprocessor
+++ b/preprocessor
@@ -1 +1 @@
-Subproject commit 6d9fc367d034be8db19857995f413fd8a740ebba
+Subproject commit 378d00fc3a50a11c3dce615ea337ad55989c938e
diff --git a/tests/analytic_derivatives/BrockMirman_PertParamsDerivs.mod b/tests/analytic_derivatives/BrockMirman_PertParamsDerivs.mod
index bcfc7735f815749b241de53ba8b8c1b73db0a798..0c5ac994dd4fa916bbfee6beb3a267583e6b4cab 100644
--- a/tests/analytic_derivatives/BrockMirman_PertParamsDerivs.mod
+++ b/tests/analytic_derivatives/BrockMirman_PertParamsDerivs.mod
@@ -418,7 +418,7 @@ for jj = 1:2
if jj==1
strparamset = 'PRIOR';
nSYM = nSYMprior;
- xparam_prior = set_prior(estim_params_,M_,options_);
+ [xparam_prior, estim_params_]= set_prior(estim_params_,M_,options_);
M_ = set_all_parameters(xparam_prior,estim_params_,M_);
elseif jj==2
strparamset = 'CALIBRATION';
diff --git a/tests/analytic_derivatives/burnside_3_order_PertParamsDerivs.mod b/tests/analytic_derivatives/burnside_3_order_PertParamsDerivs.mod
index c63091bf36e39a5971dbe4874d0c8cef52891958..ba7f17762250ed84b8c7f164739df5201cd1aa72 100644
--- a/tests/analytic_derivatives/burnside_3_order_PertParamsDerivs.mod
+++ b/tests/analytic_derivatives/burnside_3_order_PertParamsDerivs.mod
@@ -113,7 +113,7 @@ stoch_simul(order=@{ORDER},k_order_solver,irf=0,drop=0,periods=0,nograph);
identification(order=@{ORDER},nograph,no_identification_strength);
%make sure everything is computed at prior mean
-xparam_prior = set_prior(estim_params_,M_,options_);
+[xparam_prior, estim_params_]= set_prior(estim_params_,M_,options_);
M_ = set_all_parameters(xparam_prior,estim_params_,M_);
[oo_.dr,info,M_.params] = resol(0,M_, options_, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
diff --git a/tests/gsa/cod_ML_morris/cod_ML_morris.mod b/tests/gsa/cod_ML_morris/cod_ML_morris.mod
index 214556caafe479245d4f63ae3ac63e4141ba49f4..4687bd0eea5486eaa8ed5b59937f49f713e0497f 100644
--- a/tests/gsa/cod_ML_morris/cod_ML_morris.mod
+++ b/tests/gsa/cod_ML_morris/cod_ML_morris.mod
@@ -1,6 +1,6 @@
-var c, y_ro, psi, s, tt, pi_h, mc_ro, i_ro, pi_ro, pi_eu, pi_f, y_eu, mc_eu, i_eu, g, a, cp, bnr, gs, as ;
-varexo eps_a, eps_g, eps_cp, eps_s, eps_m, eps_as, eps_gs, eps_ms;
-parameters alfa, niu, delt, teta_h, teta_f, bet, fi, sigm, h, ro, psi_pi, tetas, rhos, psi_pis, psi_y, psi_deltay, rho_g, rho_a, rho_cp, rho_s, rho_gs, rho_as ;
+var c, y_ro, psi, s, tt, pi_h ${\pi_h}$, mc_ro, i_ro, pi_ro ${\pi_\rho}$, pi_eu, pi_f, y_eu, mc_eu, i_eu, g, a, cp, bnr, gs, as ;
+varexo eps_a, eps_g, eps_cp, eps_s, eps_m $\varepsilon_m$, eps_as, eps_gs, eps_ms;
+parameters alfa, niu, delt, teta_h, teta_f, bet, fi, sigm, h, ro, psi_pi $\psi_\pi$, tetas, rhos, psi_pis, psi_y, psi_deltay, rho_g, rho_a, rho_cp, rho_s, rho_gs, rho_as ;
alfa = 0.425 ;
bet = 0.99 ;
@@ -121,7 +121,7 @@ stderr eps_as, 1, 0.01, 10 ;
end;
varobs y_ro, pi_ro, i_ro, s, y_eu, pi_eu, i_eu, tt ;
-
+options_.TeX =1;
dynare_sensitivity (identification=1, nsam = 500, lik_only = 1, morris=2) ;
stoch_simul(order=2,irf=20) y_ro pi_ro i_ro s ;
diff --git a/tests/identification/forward_looking/forward_looking_empty_ghx.mod b/tests/identification/forward_looking/forward_looking_empty_ghx.mod
index 2959738c7f805aa184578b8dcd5090ed99150d31..dd22e278eaee30490cdf8b47428c8766932727a8 100755
--- a/tests/identification/forward_looking/forward_looking_empty_ghx.mod
+++ b/tests/identification/forward_looking/forward_looking_empty_ghx.mod
@@ -26,7 +26,7 @@ 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
-
+stoch_simul(order=1,noprint,irf=0,nomoments);
if max(max(abs(TRUE_SOLUTION1 - oo_.dr.ghu))) > 1e-15
error('Something wrong with perturbation');
end
diff --git a/tests/identification/kim/kim2.mod b/tests/identification/kim/kim2.mod
index bc8027d1287c3be8fa9434ce27fa949971c78c8e..54d2ae180c727b31ec6b287a7c77dec8a9197a86 100644
--- a/tests/identification/kim/kim2.mod
+++ b/tests/identification/kim/kim2.mod
@@ -125,4 +125,5 @@ end
% Integration test if identification works without priors
estim_params_=[];
+dumpy=0;
identification(advanced=1,max_dim_cova_group=3);