Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision
  • 4.3
  • 4.4
  • 4.5
  • 4.6
  • 5.x
  • 6.x
  • asm
  • aux_func
  • clang+openmp
  • dates-and-dseries-improvements
  • declare_vars_in_model_block
  • dmm
  • dragonfly
  • dynare_minreal
  • eigen
  • error_msg_undeclared_model_vars
  • estim_params
  • exo_steady_state
  • gpm-optimal-policy
  • julia
  • madysson
  • master
  • mex-GetPowerDeriv
  • penalty
  • separateM_
  • slice
  • sphinx-doc-experimental
  • static_aux_vars
  • time-varying-information-set
  • various_fixes
  • 3.062
  • 3.063
  • 4.0.0
  • 4.0.1
  • 4.0.2
  • 4.0.3
  • 4.0.4
  • 4.1-alpha1
  • 4.1-alpha2
  • 4.1.0
  • 4.1.1
  • 4.1.2
  • 4.1.3
  • 4.2.0
  • 4.2.1
  • 4.2.2
  • 4.2.3
  • 4.2.4
  • 4.2.5
  • 4.3.0
  • 4.3.1
  • 4.3.2
  • 4.3.3
  • 4.4-beta1
  • 4.4.0
  • 4.4.1
  • 4.4.2
  • 4.4.3
  • 4.5.0
  • 4.5.1
  • 4.5.2
  • 4.5.3
  • 4.5.4
  • 4.5.5
  • 4.5.6
  • 4.5.7
  • 4.6-beta1
  • 4.6.0
  • 4.6.0-rc1
  • 4.6.0-rc2
  • 4.6.1
  • 4.6.2
  • 4.6.3
  • 4.6.4
  • 4.7-beta1
  • 4.7-beta2
  • 4.7-beta3
  • 5.0
  • 5.0-rc1
  • 5.1
  • 5.2
  • 5.3
  • 5.4
  • 5.5
  • 6-beta1
  • 6-beta2
  • 6.0
  • 6.1
  • 6.2
  • 6.3
  • 6.4
91 results

Target

Select target project
  • giovanma/dynare
  • giorgiomas/dynare
  • Vermandel/dynare
  • Dynare/dynare
  • normann/dynare
  • MichelJuillard/dynare
  • wmutschl/dynare
  • FerhatMihoubi/dynare
  • sebastien/dynare
  • lnsongxf/dynare
  • rattoma/dynare
  • CIMERS/dynare
  • FredericKarame/dynare
  • SumuduK/dynare
  • MinjeJeon/dynare
  • camilomrch/dynare
  • DoraK/dynare
  • avtishin/dynare
  • selma/dynare
  • claudio_olguin/dynare
  • jeffjiang07/dynare
  • EthanSystem/dynare
  • stepan-a/dynare
  • wjgatt/dynare
  • JohannesPfeifer/dynare
  • gboehl/dynare
  • ebenetce/dynare
  • chskcau/dynare-doc-fixes
28 results
Select Git revision
  • bayesian_irf_matching
  • cet-octave-one-more-time
  • docker-meson
  • identificationME
  • irf-matching-example
  • irf_match_perfect_foresight
  • main
  • mom_fname
  • new-desktop-matlab
  • optimizers
  • pskf
  • 3.062
  • 3.063
  • 4.0.0
  • 4.0.1
  • 4.0.2
  • 4.0.3
  • 4.0.4
  • 4.1-alpha1
  • 4.1-alpha2
  • 4.1.0
  • 4.1.1
  • 4.1.2
  • 4.1.3
  • 4.2.0
  • 4.2.1
  • 4.2.2
  • 4.2.3
  • 4.2.4
  • 4.2.5
  • 4.3.0
  • 4.3.1
  • 4.3.2
  • 4.3.3
  • 4.4-beta1
  • 4.4.0
  • 4.4.1
  • 4.4.2
  • 4.4.3
  • 4.5.0
  • 4.5.1
  • 4.5.2
  • 4.5.3
  • 4.5.4
  • 4.5.5
  • 4.5.6
46 results
Show changes
Show whitespace changes
Inline Side-by-side
Showing
with 1419 additions and 485 deletions
matlab/identification_numerical_objective.m → matlab/+identification/numerical_objective.m
View file @ 4a32acc2
function out = identification_numerical_objective(params, outputflag, estim_params, M, oo, options, indpmodel, indpstderr, indpcorr, indvar, useautocorr, nlags, grid_nbr)
%function out = identification_numerical_objective(params, outputflag, estim_params, M, oo, options, indpmodel, indpstderr, indpcorr, indvar, useautocorr, nlags, grid_nbr)
function out = numerical_objective(params, outputflag, estim_params_, M_, options_, indpmodel, indpstderr, indvar, useautocorr, nlags, grid_nbr, dr, steady_state, exo_steady_state, exo_det_steady_state)
% out = numerical_objective(params, outputflag, estim_params_, M_, options_, indpmodel, indpstderr, indvar, useautocorr, nlags, grid_nbr, dr, steady_state, exo_steady_state, exo_det_steady_state)
% -------------------------------------------------------------------------
% Objective function to compute numerically the Jacobians used for identification analysis
% Previously this function was called thet2tau.m
......@@ -8,19 +8,21 @@ function out = identification_numerical_objective(params, outputflag, estim_para
% params: [vector] parameter values at which to evaluate objective function
% stderr parameters come first, corr parameters second, model parameters third
% outputflag: [integer] flag which objective to compute (see below)
% estim_params: [structure] storing the estimation information
% M: [structure] storing the model information
% oo: [structure] storing the reduced form solution results
% options: [structure] storing the options
% estim_params_: [structure] storing the estimation information
% M_: [structure] storing the model information
% options_: [structure] storing the options
% indpmodel: [vector] Index of model parameters
% indpstderr: [vector] Index of stderr parameters
% indpcorr: [matrix] Index of corr parameters
% indvar: [vector] Index of selected or observed variables
% dr [structure] Reduced form model.
% endo_steady_state [vector] steady state value for endogenous variables
% exo_steady_state [vector] steady state value for exogenous variables
% exo_det_steady_state [vector] steady state value for exogenous deterministic variables
% -------------------------------------------------------------------------
% OUTPUTS
% out: dependent on outputflag
% * 0: out = [Yss; vec(A); vec(B); dyn_vech(Sig_e)]; of indvar variables only, in DR order. This is needed to compute dTAU and Komunjer and Ng's D.
% Note that Jacobian of Om is computed in get_identification_Jacobians.m (previously getJJ.m) or get_first_order_solution_params_deriv.m (previously getH.m) from Jacobian of B and Sigma_e, because this is more efficient due to some testing with analytical derivatives from An and Schorfheide model
% Note that Jacobian of Om is computed in identification.get_jacobians.m (previously getJJ.m) or get_first_order_solution_params_deriv.m (previously getH.m) from Jacobian of B and Sigma_e, because this is more efficient due to some testing with analytical derivatives from An and Schorfheide model
% * 1: out = [vech(cov(Y_t,Y_t)); vec(cov(Y_t,Y_{t-1}); ...; vec(cov(Y_t,Y_{t-nlags})] of indvar variables, in DR order. This is needed to compute Iskrev's J.
% * 2: out = vec(spectral density) with dimension [var_nbr^2*grid_nbr,1] Spectral density of indvar variables evaluated at (grid_nbr/2+1) discretized points in the interval [0;pi]. This is needed for Qu and Tkachenko's G.
% * -1: out = g1(:); of all variables, in DR order. This is needed to compute dLRE.
......@@ -30,10 +32,10 @@ function out = identification_numerical_objective(params, outputflag, estim_para
% Jacobian of the dynamic model equations, and Y_t selected variables
% -------------------------------------------------------------------------
% This function is called by
% * get_identification_jacobians.m (previously getJJ.m)
% * identification.get_jacobians.m (previously getJJ.m)
% -------------------------------------------------------------------------
% This function calls
% * [M.fname,'.dynamic']
% * [M_.fname,'.dynamic']
% * dyn_vech
% * resol
% * vec
......@@ -59,26 +61,26 @@ function out = identification_numerical_objective(params, outputflag, estim_para
%% Update stderr, corr and model parameters
%note that if no estimated_params_block is given, then all stderr and model parameters are selected but no corr parameters
if length(params) > length(indpmodel)
if isempty(indpstderr)==0 && isempty(estim_params.var_exo) %if there are stderr parameters but no estimated_params_block
if isempty(indpstderr)==0 && isempty(estim_params_.var_exo) %if there are stderr parameters but no estimated_params_block
%provide temporary necessary information for stderr parameters
estim_params.nvx = length(indpstderr);
estim_params.var_exo = indpstderr';
estim_params_.nvx = length(indpstderr);
estim_params_.var_exo = indpstderr';
end
if isempty(indpmodel)==0 && isempty(estim_params.param_vals) %if there are model parameters but no estimated_params_block
if isempty(indpmodel)==0 && isempty(estim_params_.param_vals) %if there are model parameters but no estimated_params_block
%provide temporary necessary information for model parameters
estim_params.np = length(indpmodel);
estim_params.param_vals = indpmodel';
estim_params_.np = length(indpmodel);
estim_params_.param_vals = indpmodel';
end
M = set_all_parameters(params,estim_params,M); %this function can only be used if there is some information in estim_params
M_ = set_all_parameters(params,estim_params_,M_); %this function can only be used if there is some information in estim_params_
else
%if there are only model parameters, we don't need to use set_all_parameters
M.params(indpmodel) = params;
M_.params(indpmodel) = params;
end
%% compute Kalman transition matrices and steady state with updated parameters
[oo.dr,info,M.params] = compute_decision_rules(M,options,oo.dr, oo.steady_state, oo.exo_steady_state, oo.exo_det_steady_state);
options = rmfield(options,'options_ident');
pruned = pruned_state_space_system(M, options, oo.dr, indvar, nlags, useautocorr, 0);
[dr,~,M_.params] = compute_decision_rules(M_,options_,dr, steady_state, exo_steady_state, exo_det_steady_state);
options_ = rmfield(options_,'options_ident');
pruned = pruned_SS.pruned_state_space_system(M_, options_, dr, indvar, nlags, useautocorr, 0);
%% out = [vech(cov(Y_t,Y_t)); vec(cov(Y_t,Y_{t-1}); ...; vec(cov(Y_t,Y_{t-nlags})] of indvar variables, in DR order. This is Iskrev (2010)'s J matrix.
if outputflag == 1
......
matlab/plot_identification.m → matlab/+identification/plot.m
View file @ 4a32acc2
function plot_identification(params, idemoments, idehess, idemodel, idelre, advanced, tittxt, name, IdentifDirectoryName, tit_TeX, name_tex)
% function plot_identification(params,idemoments,idehess,idemodel, idelre, advanced, tittxt, name, IdentifDirectoryName)
function plot(M_, params, idemoments, idehess, idemodel, idelre, advanced, tittxt, name, IdentifDirectoryName, fname, options_, estim_params_, bayestopt_, tit_TeX, name_tex)
% plot(M_, 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
......@@ -9,17 +10,22 @@ function plot_identification(params, idemoments, idehess, idemodel, idelre, adva
% o idelre [structure] identification results for the LRE model
% o advanced [integer] flag for advanced identification checks
% o tittxt [char] name of the results to plot
% o name [char] list of names
% o name [char] list of parameter names
% o IdentifDirectoryName [char] directory name
% o fname [char] file name
% o options_ [structure] structure describing the current options
% o estim_params_ [structure] characterizing parameters to be estimated
% o bayestopt_ [structure] describing the priors
% o tittxt [char] TeX-name of the results to plot
% o name_tex [char] TeX-names of the parameters
%
% OUTPUTS
% None
%
% SPECIAL REQUIREMENTS
% None
% Copyright © 2008-2023 Dynare Team
% Copyright © 2008-2025 Dynare Team
%
% This file is part of Dynare.
%
......@@ -36,14 +42,12 @@ 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/>.
global M_ options_
if nargin <10 || isempty(tit_TeX)
if nargin <14 || isempty(tit_TeX)
tit_TeX=tittxt;
end
if nargin <11
name_TeX=name;
if nargin <15
name_tex=name;
end
[SampleSize, nparam]=size(params);
......@@ -54,20 +58,18 @@ si_dLREnorm = idelre.si_dDYNAMICnorm;
tittxt1=regexprep(tittxt, ' ', '_');
tittxt1=strrep(tittxt1, '.', '');
if SampleSize == 1
si_dMOMENTS = idemoments.si_dMOMENTS;
hh_fig = dyn_figure(options_.nodisplay,'Name',[tittxt, ' - Identification using info from observables']);
subplot(211)
mmm = (idehess.ide_strength_dMOMENTS);
[ss, is] = sort(mmm);
if ~all(isnan(idehess.ide_strength_dMOMENTS_prior)) ...
&& ~(nparam == 1 && ~isoctave && matlab_ver_less_than('9.7')) % MATLAB < R2019b does not accept bar(1, [2 3])
[~, is] = sort(mmm);
if ~all(isnan(idehess.ide_strength_dMOMENTS_prior))
bar(1:nparam,log([idehess.ide_strength_dMOMENTS(:,is)' idehess.ide_strength_dMOMENTS_prior(:,is)']))
else
bar(1:nparam,log([idehess.ide_strength_dMOMENTS(:,is)' ]))
bar(1:nparam,log(idehess.ide_strength_dMOMENTS(:,is)'))
end
hold on
plot((1:length(idehess.ide_strength_dMOMENTS(:,is)))-0.15,log([idehess.ide_strength_dMOMENTS(:,is)']),'o','MarkerSize',7,'MarkerFaceColor',[0 0 0],'MarkerEdgeColor','none')
plot((1:length(idehess.ide_strength_dMOMENTS_prior(:,is)))+0.15,log([idehess.ide_strength_dMOMENTS_prior(:,is)']),'o','MarkerSize',7,'MarkerFaceColor',[0 0 0],'MarkerEdgeColor','none')
plot((1:length(idehess.ide_strength_dMOMENTS(:,is)))-0.15,log(idehess.ide_strength_dMOMENTS(:,is)'),'o','MarkerSize',7,'MarkerFaceColor',[0 0 0],'MarkerEdgeColor','none')
plot((1:length(idehess.ide_strength_dMOMENTS_prior(:,is)))+0.15,log(idehess.ide_strength_dMOMENTS_prior(:,is)'),'o','MarkerSize',7,'MarkerFaceColor',[0 0 0],'MarkerEdgeColor','none')
if any(isinf(log(idehess.ide_strength_dMOMENTS(idehess.identified_parameter_indices))))
%-Inf, i.e. 0 strength
inf_indices=find(isinf(log(idehess.ide_strength_dMOMENTS(idehess.identified_parameter_indices))) & log(idehess.ide_strength_dMOMENTS(idehess.identified_parameter_indices))<0);
......@@ -92,8 +94,12 @@ if SampleSize == 1
set(gca,'xticklabel','')
dy = get(gca,'ylim');
for ip=1:nparam
if options_.TeX
text(ip,dy(1),name_tex{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','latex')
else
text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
end
if ~all(isnan(idehess.ide_strength_dMOMENTS_prior))
legend('relative to param value','relative to prior std','Location','Best')
else
......@@ -106,15 +112,14 @@ if SampleSize == 1
end
subplot(212)
if ~all(isnan(idehess.deltaM_prior)) ...
&& ~(nparam == 1 && ~isoctave && matlab_ver_less_than('9.7')) % MATLAB < R2019b does not accept bar(1, [2 3])
if ~all(isnan(idehess.deltaM_prior))
bar(1:nparam, log([idehess.deltaM(is) idehess.deltaM_prior(is)]))
else
bar(1:nparam, log([idehess.deltaM(is)]))
end
hold on
plot((1:length(idehess.deltaM(is)))-0.15,log([idehess.deltaM(is)']),'o','MarkerSize',7,'MarkerFaceColor',[0 0 0],'MarkerEdgeColor','none')
plot((1:length(idehess.deltaM_prior(is)))+0.15,log([idehess.deltaM_prior(is)']),'o','MarkerSize',7,'MarkerFaceColor',[0 0 0],'MarkerEdgeColor','none')
plot((1:length(idehess.deltaM(is)))-0.15,log(idehess.deltaM(is)'),'o','MarkerSize',7,'MarkerFaceColor',[0 0 0],'MarkerEdgeColor','none')
plot((1:length(idehess.deltaM_prior(is)))+0.15,log(idehess.deltaM_prior(is)'),'o','MarkerSize',7,'MarkerFaceColor',[0 0 0],'MarkerEdgeColor','none')
inf_pos=find(isinf(log(idehess.deltaM)));
if ~isempty(inf_pos)
inf_indices=~ismember(inf_pos,idehess.sensitivity_zero_pos);
......@@ -131,8 +136,12 @@ if SampleSize == 1
set(gca,'xticklabel','')
dy = get(gca,'ylim');
for ip=1:nparam
if options_.TeX
text(ip,dy(1),name_tex{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','latex')
else
text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
end
if ~all(isnan(idehess.deltaM_prior))
legend('relative to param value','relative to prior std','Location','Best')
else
......@@ -144,19 +153,19 @@ if SampleSize == 1
title('Sensitivity component with moments Information matrix (log-scale)')
end
if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
fidTeX = fopen([IdentifDirectoryName '/' M_.fname '_ident_strength_' tittxt1,'.tex'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by plot_identification.m (Dynare).\n');
fidTeX = fopen([IdentifDirectoryName '/' fname '_ident_strength_' tittxt1,'.tex'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by identification.plot.m (Dynare).\n');
fprintf(fidTeX,['%% ' datestr(now,0) '\n\n']);
fprintf(fidTeX,'\\begin{figure}[H]\n');
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[width=0.8\\textwidth]{%s_ident_strength_%s}\n',[IdentifDirectoryName '/' M_.fname],tittxt1);
fprintf(fidTeX,'\\includegraphics[width=0.8\\textwidth]{%s_ident_strength_%s}\n',[IdentifDirectoryName '/' fname],tittxt1);
fprintf(fidTeX,'\\caption{%s - Identification using info from observables.}',tit_TeX);
fprintf(fidTeX,'\\label{Fig:ident:%s}\n',deblank(tittxt));
fprintf(fidTeX,'\\end{figure}\n\n');
fprintf(fidTeX,'%% End Of TeX file. \n');
fclose(fidTeX);
end
dyn_saveas(hh_fig,[IdentifDirectoryName '/' M_.fname '_ident_strength_' tittxt1],options_.nodisplay,options_.graph_format);
dyn_saveas(hh_fig,[IdentifDirectoryName '/' fname '_ident_strength_' tittxt1],options_.nodisplay,options_.graph_format);
if advanced
if ~options_.nodisplay
......@@ -181,18 +190,22 @@ if SampleSize == 1
set(gca,'xticklabel','')
dy = get(gca,'ylim');
for ip=1:nparam
if options_.TeX
text(ip,dy(1),name_tex{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','latex')
else
text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
end
legend('Moments','Model','LRE model','Location','Best')
title('Sensitivity bars using derivatives (log-scale)')
dyn_saveas(hh_fig,[IdentifDirectoryName '/' M_.fname '_sensitivity_' tittxt1 ],options_.nodisplay,options_.graph_format);
dyn_saveas(hh_fig,[IdentifDirectoryName '/' fname '_sensitivity_' tittxt1 ],options_.nodisplay,options_.graph_format);
if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
fidTeX = fopen([IdentifDirectoryName '/' M_.fname '_sensitivity_' tittxt1,'.tex'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by plot_identification.m (Dynare).\n');
fidTeX = fopen([IdentifDirectoryName '/' fname '_sensitivity_' tittxt1,'.tex'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by identification.plot.m (Dynare).\n');
fprintf(fidTeX,['%% ' datestr(now,0) '\n\n']);
fprintf(fidTeX,'\\begin{figure}[H]\n');
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[width=0.8\\textwidth]{%s_sensitivity_%s}\n',[IdentifDirectoryName '/' M_.fname],tittxt1);
fprintf(fidTeX,'\\includegraphics[width=0.8\\textwidth]{%s_sensitivity_%s}\n',[IdentifDirectoryName '/' fname],tittxt1);
fprintf(fidTeX,'\\caption{%s - Sensitivity plot.}',tit_TeX);
fprintf(fidTeX,'\\label{Fig:sensitivity:%s}\n',deblank(tittxt));
fprintf(fidTeX,'\\end{figure}\n\n');
......@@ -203,30 +216,35 @@ if SampleSize == 1
% identificaton patterns
for j=1:size(idemoments.cosndMOMENTS,2)
pax=NaN(nparam,nparam);
% fprintf('\n')
% disp(['Collinearity patterns with ', int2str(j) ,' parameter(s)'])
% fprintf('%-15s [%-*s] %10s\n','Parameter',(15+1)*j,' Expl. params ','cosn')
for i=1:nparam
namx='';
for in=1:j
dumpindx = idemoments.pars{i,j}(in);
if isnan(dumpindx)
namx=[namx ' ' sprintf('%-15s','--')];
else
if options_.TeX
namx=[namx ' ' sprintf('%-15s',name_tex{dumpindx})];
else
namx=[namx ' ' sprintf('%-15s',name{dumpindx})];
end
pax(i,dumpindx)=idemoments.cosndMOMENTS(i,j);
end
end
% fprintf('%-15s [%s] %10.3f\n',name{i},namx,idemoments.cosndMOMENTS(i,j))
end
hh_fig = dyn_figure(options_.nodisplay,'Name',[tittxt,' - Collinearity patterns with ', int2str(j) ,' parameter(s)']);
imagesc(pax,[0 1]);
set(gca,'xticklabel','')
set(gca,'yticklabel','')
for ip=1:nparam
if options_.TeX
text(ip,(0.5),name_tex{ip},'rotation',90,'HorizontalAlignment','left','interpreter','latex')
text(0.5,ip,name_tex{ip},'rotation',0,'HorizontalAlignment','right','interpreter','latex')
else
text(ip,(0.5),name{ip},'rotation',90,'HorizontalAlignment','left','interpreter','none')
text(0.5,ip,name{ip},'rotation',0,'HorizontalAlignment','right','interpreter','none')
end
end
colorbar;
colormap('jet');
ax=colormap;
......@@ -239,14 +257,14 @@ if SampleSize == 1
set(gca,'xgrid','on')
set(gca,'ygrid','on')
xlabel([tittxt,' - Collinearity patterns with ', int2str(j) ,' parameter(s)'],'interpreter','none')
dyn_saveas(hh_fig,[ IdentifDirectoryName '/' M_.fname '_ident_collinearity_' tittxt1 '_' int2str(j) ],options_.nodisplay,options_.graph_format);
dyn_saveas(hh_fig,[ IdentifDirectoryName '/' fname '_ident_collinearity_' tittxt1 '_' int2str(j) ],options_.nodisplay,options_.graph_format);
if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
fidTeX = fopen([ IdentifDirectoryName '/' M_.fname '_ident_collinearity_' tittxt1 '_' int2str(j),'.tex'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by plot_identification.m (Dynare).\n');
fidTeX = fopen([ IdentifDirectoryName '/' fname '_ident_collinearity_' tittxt1 '_' int2str(j),'.tex'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by identification.plot.m (Dynare).\n');
fprintf(fidTeX,['%% ' datestr(now,0) '\n\n']);
fprintf(fidTeX,'\\begin{figure}[H]\n');
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[width=0.8\\textwidth]{%s_ident_collinearity_%s_%u}\n',[IdentifDirectoryName '/' M_.fname],tittxt1,j);
fprintf(fidTeX,'\\includegraphics[width=0.8\\textwidth]{%s_ident_collinearity_%s_%u}\n',[IdentifDirectoryName '/' fname],tittxt1,j);
fprintf(fidTeX,'\\caption{%s - Collinearity patterns with %u parameter(s).}',tit_TeX,j);
fprintf(fidTeX,'\\label{Fig:collinearity:%s:%u_pars}\n',deblank(tittxt),j);
fprintf(fidTeX,'\\end{figure}\n\n');
......@@ -255,7 +273,7 @@ if SampleSize == 1
end
end
skipline()
[U,S,V]=svd(idehess.AHess,0);
[~,S,V]=svd(idehess.AHess,0);
S=diag(S);
if idehess.flag_score
if nparam<5
......@@ -268,8 +286,6 @@ if SampleSize == 1
tex_tit_2=[tittxt,' - Identification patterns (Information matrix): HIGHEST SV'];
end
else
% S = idemoments.S;
% V = idemoments.V;
if nparam<5
f1 = dyn_figure(options_.nodisplay,'Name',[tittxt,' - Identification patterns (moments Information matrix)']);
tex_tit_1=[tittxt,' - Identification patterns (moments Information matrix)'];
......@@ -299,19 +315,23 @@ if SampleSize == 1
set(gca,'xticklabel','')
if j==4 || j==nparam || j==8
for ip=1:nparam
if options_.TeX
text(ip,-0.02,name_tex{ip},'rotation',90,'HorizontalAlignment','right','interpreter','latex')
else
text(ip,-0.02,name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
end
end
title(['Singular value ',num2str(Stit)])
end
dyn_saveas(f1,[ IdentifDirectoryName '/' M_.fname '_ident_pattern_' tittxt1 '_1' ],options_.nodisplay,options_.graph_format);
dyn_saveas(f1,[ IdentifDirectoryName '/' fname '_ident_pattern_' tittxt1 '_1' ],options_.nodisplay,options_.graph_format);
if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
fidTeX = fopen([ IdentifDirectoryName '/' M_.fname '_ident_pattern_' tittxt1 '_1','.tex'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by plot_identification.m (Dynare).\n');
fidTeX = fopen([ IdentifDirectoryName '/' fname '_ident_pattern_' tittxt1 '_1','.tex'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by identification.plot.m (Dynare).\n');
fprintf(fidTeX,['%% ' datestr(now,0) '\n\n']);
fprintf(fidTeX,'\\begin{figure}[H]\n');
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[width=0.8\\textwidth]{%s_ident_pattern_%s_1}\n',[IdentifDirectoryName '/' M_.fname],tittxt1);
fprintf(fidTeX,'\\includegraphics[width=0.8\\textwidth]{%s_ident_pattern_%s_1}\n',[IdentifDirectoryName '/' fname],tittxt1);
fprintf(fidTeX,'\\caption{%s.}',tex_tit_1);
fprintf(fidTeX,'\\label{Fig:ident_pattern:%s:1}\n',tittxt1);
fprintf(fidTeX,'\\end{figure}\n\n');
......@@ -319,14 +339,14 @@ if SampleSize == 1
fclose(fidTeX);
end
if nparam>4
dyn_saveas(f2,[ IdentifDirectoryName '/' M_.fname '_ident_pattern_' tittxt1 '_2' ],options_.nodisplay,options_.graph_format);
dyn_saveas(f2,[ IdentifDirectoryName '/' fname '_ident_pattern_' tittxt1 '_2' ],options_.nodisplay,options_.graph_format);
if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
fidTeX = fopen([ IdentifDirectoryName '/' M_.fname '_ident_pattern_' tittxt1 '_2.tex'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by plot_identification.m (Dynare).\n');
fidTeX = fopen([ IdentifDirectoryName '/' fname '_ident_pattern_' tittxt1 '_2.tex'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by identification.plot.m (Dynare).\n');
fprintf(fidTeX,['%% ' datestr(now,0) '\n\n']);
fprintf(fidTeX,'\\begin{figure}[H]\n');
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[width=0.8\\textwidth]{%s_ident_pattern_%s_2}\n',[IdentifDirectoryName '/' M_.fname],tittxt1);
fprintf(fidTeX,'\\includegraphics[width=0.8\\textwidth]{%s_ident_pattern_%s_2}\n',[IdentifDirectoryName '/' fname],tittxt1);
fprintf(fidTeX,'\\caption{%s.}',tex_tit_2);
fprintf(fidTeX,'\\label{Fig:ident_pattern:%s:2}\n',tittxt1);
fprintf(fidTeX,'\\end{figure}\n\n');
......@@ -337,10 +357,10 @@ if SampleSize == 1
end
else
hh_fig = dyn_figure(options_.nodisplay,'Name',['MC sensitivities']);
hh_fig = dyn_figure(options_.nodisplay,'Name','MC sensitivities');
subplot(211)
mmm = (idehess.ide_strength_dMOMENTS);
[ss, is] = sort(mmm);
[~, is] = sort(mmm);
mmm = mean(si_dMOMENTSnorm)';
mmm = mmm./max(mmm);
if advanced
......@@ -357,20 +377,24 @@ else
set(gca,'xticklabel','')
dy = get(gca,'ylim');
for ip=1:nparam
if options_.TeX
text(ip,dy(1),name_tex{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','latex')
else
text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
end
if advanced
legend('Moments','Model','LRE model','Location','Best')
end
title('MC mean of sensitivity measures')
dyn_saveas(hh_fig,[ IdentifDirectoryName '/' M_.fname '_MC_sensitivity' ],options_.nodisplay,options_.graph_format);
dyn_saveas(hh_fig,[ IdentifDirectoryName '/' fname '_MC_sensitivity' ],options_.nodisplay,options_.graph_format);
if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
fidTeX = fopen([ IdentifDirectoryName '/' M_.fname '_MC_sensitivity.tex'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by plot_identification.m (Dynare).\n');
fidTeX = fopen([ IdentifDirectoryName '/' fname '_MC_sensitivity.tex'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by identification.plot.m (Dynare).\n');
fprintf(fidTeX,['%% ' datestr(now,0) '\n\n']);
fprintf(fidTeX,'\\begin{figure}[H]\n');
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[width=0.8\\textwidth]{%s_MC_sensitivity}\n',[IdentifDirectoryName '/' M_.fname]);
fprintf(fidTeX,'\\includegraphics[width=0.8\\textwidth]{%s_MC_sensitivity}\n',[IdentifDirectoryName '/' fname]);
fprintf(fidTeX,'\\caption{MC mean of sensitivity measures}');
fprintf(fidTeX,'\\label{Fig:_MC_sensitivity}\n');
fprintf(fidTeX,'\\end{figure}\n\n');
......@@ -379,57 +403,62 @@ else
end
if advanced
if ~options_.nodisplay,
if ~options_.nodisplay
skipline()
disp('Displaying advanced diagnostics')
end
% options_.nograph=1;
hh_fig = dyn_figure(options_.nodisplay,'Name','MC Condition Number');
subplot(221)
if isoctave
hist(log10(idemodel.cond))
else
histogram(log10(idemodel.cond))
end
title('log10 of Condition number in the model')
subplot(222)
if isoctave
hist(log10(idemoments.cond))
else
histogram(log10(idemoments.cond))
end
title('log10 of Condition number in the moments')
subplot(223)
if isoctave
hist(log10(idelre.cond))
else
histogram(log10(idelre.cond))
end
title('log10 of Condition number in the LRE model')
dyn_saveas(hh_fig,[IdentifDirectoryName '/' M_.fname '_ident_COND' ],options_.nodisplay,options_.graph_format);
dyn_saveas(hh_fig,[IdentifDirectoryName '/' fname '_ident_COND' ],options_.nodisplay,options_.graph_format);
options_mcf.pvalue_ks = 0.1;
options_mcf.pvalue_corr = 0.001;
options_mcf.alpha2 = 0;
options_mcf.param_names = name;
options_mcf.param_names_tex = name_tex;
options_mcf.fname_ = M_.fname;
options_mcf.fname_ = fname;
options_mcf.OutputDirectoryName = IdentifDirectoryName;
options_mcf.beha_title = 'LOW condition nbr';
options_mcf.nobeha_title = 'HIGH condition nbr';
if options_.TeX
options_mcf.beha_title_latex = 'LOW condition nbr';
options_mcf.nobeha_title_latex = 'HIGH condition nbr';
end
options_mcf.amcf_name = 'MC_HighestCondNumberLRE';
options_mcf.amcf_title = 'MC Highest Condition Number LRE Model';
options_mcf.title = 'MC Highest Condition Number LRE Model';
ncut=floor(SampleSize/10*9);
[dum,is]=sort(idelre.cond);
mcf_analysis(params, is(1:ncut), is(ncut+1:end), options_mcf, options_);
[~,is]=sort(idelre.cond);
gsa.monte_carlo_filtering_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';
[dum,is]=sort(idemodel.cond);
mcf_analysis(params, is(1:ncut), is(ncut+1:end), options_mcf, options_);
[~,is]=sort(idemodel.cond);
gsa.monte_carlo_filtering_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';
[dum,is]=sort(idemoments.cond);
mcf_analysis(params, is(1:ncut), is(ncut+1:end), options_mcf, options_);
% [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);
% % inonbeh=find(idemoments.Mco(j,:)>=0.9);
% % if ~isempty(ibeh) && ~isempty(inonbeh)
% % [proba, dproba] = stab_map_1(params, ibeh, inonbeh, ['HighestMultiCollinearity_',name{j}], 1, [], IdentifDirectoryName);
% % end
% [~,is]=sort(idemoments.Mco(:,j));
% [proba, dproba] = stab_map_1(params, is(1:ncut), is(ncut+1:end), ['MC_HighestMultiCollinearity_',name{j}], 1, [], IdentifDirectoryName, 0.15);
% end
[~,is]=sort(idemoments.cond);
gsa.monte_carlo_filtering_analysis(params, is(1:ncut), is(ncut+1:end), options_mcf, M_, options_, bayestopt_, estim_params_);
if nparam<5
f1 = dyn_figure(options_.nodisplay,'Name',[tittxt,' - MC Identification patterns (moments): HIGHEST SV']);
......@@ -459,8 +488,10 @@ else
SSS = idemoments.S(:,jj);
end
subplot(nsubplo,1,jj)
post_median=NaN(1,nparam);
hpd_interval=NaN(nparam,2);
for i=1:nparam
[post_mean, post_median(:,i), post_var, hpd_interval(i,:), post_deciles] = posterior_moments(VVV(:,i),0,0.9);
[~, post_median(:,i), ~, hpd_interval(i,:)] = posterior_moments(VVV(:,i),0.9);
end
bar(post_median)
hold on, plot(hpd_interval,'--*r'),
......@@ -469,19 +500,23 @@ else
set(gca,'xticklabel','')
if j==4 || j==nparam || j==8
for ip=1:nparam
if options_.TeX
text(ip,-0.02,name_tex{ip},'rotation',90,'HorizontalAlignment','right','interpreter','latex')
else
text(ip,-0.02,name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
end
end
title(['MEAN Singular value ',num2str(Stit)])
end
dyn_saveas(f1,[IdentifDirectoryName '/' M_.fname '_MC_ident_pattern_1' ],options_.nodisplay,options_.graph_format);
dyn_saveas(f1,[IdentifDirectoryName '/' fname '_MC_ident_pattern_1' ],options_.nodisplay,options_.graph_format);
if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
fidTeX = fopen([IdentifDirectoryName '/' M_.fname '_MC_ident_pattern_1.tex'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by plot_identification.m (Dynare).\n');
fidTeX = fopen([IdentifDirectoryName '/' fname '_MC_ident_pattern_1.tex'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by identification.plot.m (Dynare).\n');
fprintf(fidTeX,['%% ' datestr(now,0) '\n\n']);
fprintf(fidTeX,'\\begin{figure}[H]\n');
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[width=0.8\\textwidth]{%s_MC_ident_pattern_1}\n',[IdentifDirectoryName '/' M_.fname]);
fprintf(fidTeX,'\\includegraphics[width=0.8\\textwidth]{%s_MC_ident_pattern_1}\n',[IdentifDirectoryName '/' fname]);
fprintf(fidTeX,'\\caption{%s.}',tex_tit_1);
fprintf(fidTeX,'\\label{Fig:MC_ident_pattern:1}\n');
fprintf(fidTeX,'\\end{figure}\n\n');
......@@ -489,14 +524,14 @@ else
fclose(fidTeX);
end
if nparam>4
dyn_saveas(f2,[ IdentifDirectoryName '/' M_.fname '_MC_ident_pattern_2' ],options_.nodisplay,options_.graph_format);
dyn_saveas(f2,[ IdentifDirectoryName '/' fname '_MC_ident_pattern_2' ],options_.nodisplay,options_.graph_format);
if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
fidTeX = fopen([ IdentifDirectoryName '/' M_.fname '_MC_ident_pattern_2.tex'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by plot_identification.m (Dynare).\n');
fidTeX = fopen([ IdentifDirectoryName '/' fname '_MC_ident_pattern_2.tex'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by identification.plot.m (Dynare).\n');
fprintf(fidTeX,['%% ' datestr(now,0) '\n\n']);
fprintf(fidTeX,'\\begin{figure}[H]\n');
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[width=0.8\\textwidth]{%s_MC_ident_pattern_2}\n',[IdentifDirectoryName '/' M_.fname]);
fprintf(fidTeX,'\\includegraphics[width=0.8\\textwidth]{%s_MC_ident_pattern_2}\n',[IdentifDirectoryName '/' fname]);
fprintf(fidTeX,'\\caption{%s.}',tex_tit_2);
fprintf(fidTeX,'\\label{Fig:MC_ident_pattern:2}\n');
fprintf(fidTeX,'\\end{figure}\n\n');
......
matlab/dynare_identification.m → matlab/+identification/run.m
View file @ 4a32acc2
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(options_ident, pdraws0)
function [pdraws, STO_REDUCEDFORM, STO_MOMENTS, STO_DYNAMIC, STO_si_dDYNAMIC, STO_si_dREDUCEDFORM, STO_si_dMOMENTS, STO_dSPECTRUM, STO_dMINIMAL] = run(M_,oo_,options_,bayestopt_,estim_params_,options_ident, pdraws0)
% [pdraws, STO_REDUCEDFORM, STO_MOMENTS, STO_DYNAMIC, STO_si_dDYNAMIC, STO_si_dREDUCEDFORM, STO_si_dMOMENTS, STO_dSPECTRUM, STO_dMINIMAL] = run(options_ident, pdraws0)
% -------------------------------------------------------------------------
% This function is called, when the user specifies identification(...); in the mod file. It prepares all identification analysis:
% (1) set options, local and persistent variables for a new identification
......@@ -11,6 +11,11 @@ function [pdraws, STO_REDUCEDFORM, STO_MOMENTS, STO_DYNAMIC, STO_si_dDYNAMIC, ST
% to put identification in your mod file, otherwise the preprocessor won't provide all necessary objects
% =========================================================================
% 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_ident [structure] identification options
% * pdraws0 [SampleSize by totparam_nbr] optional: matrix of MC sample of model parameters
% -------------------------------------------------------------------------
......@@ -27,19 +32,19 @@ function [pdraws, STO_REDUCEDFORM, STO_MOMENTS, STO_DYNAMIC, STO_si_dDYNAMIC, ST
% -------------------------------------------------------------------------
% This function is called by
% * driver.m
% * map_ident_.m
% * gsa.map_identification.m
% -------------------------------------------------------------------------
% This function calls
% * checkpath
% * disp_identification
% * identification.display
% * dyn_waitbar
% * dyn_waitbar_close
% * get_all_parameters
% * get_posterior_parameters
% * get_the_name
% * identification_analysis
% * identification.analysis
% * isoctave
% * plot_identification
% * identification.plot
% * dprior.draw
% * set_default_option
% * set_prior
......@@ -65,26 +70,11 @@ 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_
% The TeX option crashes MATLAB R2014a run with "-nodisplay" option
% (as is done from the testsuite).
% Since we can’t directly test whether "-nodisplay" has been passed,
% we test for the "TOP_TEST_DIR" environment variable, which is set
% by the testsuite.
% Note that it was not tested whether the crash happens with more
% recent MATLAB versions, so when OLD_MATLAB_VERSION is increased,
% one should make a test before removing this workaround.
if options_.TeX && ~isoctave && matlab_ver_less_than('8.4') && ~isempty(getenv('TOP_TEST_DIR'))
warning('Disabling TeX option due to a bug in MATLAB R2014a with -nodisplay')
options_.TeX = false;
end
store_options_ = options_; % store options to restore them at the end
fname = M_.fname; %model name
dname = M_.dname; %model name
%turn warnings off, either globally or only relevant ids
orig_warning_state = warning;
if isoctave
%warning('off'),
warning('off','Octave:singular-matrix');
......@@ -106,7 +96,7 @@ end
options_ident = set_default_option(options_ident,'gsa_sample_file',0);
% 0: do not use sample file
% 1: triggers gsa prior sample
% 2: triggers gsa Monte-Carlo sample (i.e. loads a sample corresponding to pprior=0 and ppost=0 in dynare_sensitivity options)
% 2: triggers gsa Monte-Carlo sample (i.e. loads a sample corresponding to pprior=0 and ppost=0 in sensitivity.run options)
% FILENAME: use sample file in provided path
options_ident = set_default_option(options_ident,'parameter_set','prior_mean');
% 'calibration': use values in M_.params and M_.Sigma_e to update estimated stderr, corr and model parameters (get_all_parameters)
......@@ -151,7 +141,7 @@ options_ident = set_default_option(options_ident,'tol_rank','robust');
options_ident = set_default_option(options_ident,'tol_deriv',1.e-8);
% tolerance level for selecting columns of non-zero derivatives
options_ident = set_default_option(options_ident,'tol_sv',1.e-3);
% tolerance level for selecting non-zero singular values in identification_checks.m
% tolerance level for selecting non-zero singular values in identification.checks.m
options_ident = set_default_option(options_ident,'schur_vec_tol',1e-11);
% tolerance level used to find nonstationary variables in Schur decomposition of the transition matrix.
......@@ -192,7 +182,7 @@ if (isfield(options_ident,'no_identification_strength') && options_ident.no_ide
options_ident.no_identification_moments = 0;
end
%overwrite setting, as dynare_sensitivity does not make use of spectrum and minimal system
%overwrite setting, as sensitivity.run does not make use of spectrum and minimal system
if isfield(options_,'opt_gsa') && isfield(options_.opt_gsa,'identification') && options_.opt_gsa.identification == 1
options_ident.no_identification_minimal = 1;
options_ident.no_identification_spectrum = 1;
......@@ -220,7 +210,7 @@ end
% 5: i) option 2 for non-stationary elements by setting their initial variance in the forecast error matrix to 10 on the diagonal and all co-variances to 0 and
% ii) option 1 for the stationary elements
options_ident = set_default_option(options_ident,'analytic_derivation',1);
% 1: analytic derivation of gradient and hessian of likelihood in dsge_likelihood.m, only works for stationary models, i.e. kalman_algo<3
% 1: analytic derivation of gradient and Hessian of likelihood in dsge_likelihood.m, only works for stationary models, i.e. kalman_algo<3
options_ident = set_default_option(options_ident,'order',1);
% 1: first-order perturbation approximation, identification is based on linear state space system
% 2: second-order perturbation approximation, identification is based on second-order pruned state space system
......@@ -250,7 +240,7 @@ end
% check for external draws, i.e. set pdraws0 for a gsa analysis
if options_ident.gsa_sample_file
GSAFolder = checkpath('gsa',dname);
GSAFolder = CheckPath('gsa',dname);
if options_ident.gsa_sample_file==1
load([GSAFolder,filesep,fname,'_prior'],'lpmat','lpmat0','istable');
elseif options_ident.gsa_sample_file==2
......@@ -268,11 +258,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;
......@@ -307,24 +297,29 @@ options_.prior_mc = options_ident.prior_mc;
options_.schur_vec_tol = options_ident.schur_vec_tol;
options_.nomoments = 0;
options_.analytic_derivation=options_ident.analytic_derivation;
% 1: analytic derivation of gradient and hessian of likelihood in dsge_likelihood.m, only works for stationary models, i.e. kalman_algo<3
% 1: analytic derivation of gradient and Hessian of likelihood in dsge_likelihood.m, only works for stationary models, i.e. kalman_algo<3
options_ = set_default_option(options_,'datafile','');
options_.mode_compute = 0;
if strcmp('slice',options_.posterior_sampler_options.posterior_sampling_method)
if strfind(options_.posterior_sampler_options.sampling_opt,'slice_initialize_with_mode'',1')
options_.posterior_sampler_options.sampling_opt=strrep(options_.posterior_sampler_options.sampling_opt,'slice_initialize_with_mode'',1','slice_initialize_with_mode'',0') % reset option to prevent crash in check_posterior_sampler_options.m if mode_compute is set to 0, see https://git.dynare.org/Dynare/dynare/-/issues/1957
end
end
options_.plot_priors = 0;
options_.smoother = 1;
options_.options_ident = [];
[dataset_, dataset_info, xparam1, hh, M_, options_, oo_, estim_params_, bayestopt_, bounds] = dynare_estimation_init(M_.endo_names, fname, 1, M_, options_, oo_, estim_params_, bayestopt_);
[~, dataset_info, ~, ~, M_, options_, oo_, estim_params_, bayestopt_] = dynare_estimation_init(M_.endo_names, fname, 1, M_, options_, oo_, estim_params_, bayestopt_);
% set method to compute identification Jacobians (kronflag). Default:0
options_ident = set_default_option(options_ident,'analytic_derivation_mode', options_.analytic_derivation_mode); % if not set by user, inherit default global one
% 0: efficient sylvester equation method to compute analytical derivatives as in Ratto & Iskrev (2012)
% 1: kronecker products method to compute analytical derivatives as in Iskrev (2010) (only for order=1)
% -1: numerical two-sided finite difference method to compute numerical derivatives of all identification Jacobians using function identification_numerical_objective.m (previously thet2tau.m)
% -1: numerical two-sided finite difference method to compute numerical derivatives of all identification Jacobians using function identification.numerical_objective.m (previously thet2tau.m)
% -2: numerical two-sided finite difference method to compute numerically dYss, dg1, dg2, dg3, d2Yss and d2g1, the identification Jacobians are then computed analytically as with 0
if options_.discretionary_policy || options_.ramsey_policy
if options_ident.analytic_derivation_mode~=-1
fprintf('dynare_identification: discretionary_policy and ramsey_policy require analytic_derivation_mode=-1. Resetting the option.')
fprintf('identification.run: discretionary_policy and ramsey_policy require analytic_derivation_mode=-1. Resetting the option.')
options_ident.analytic_derivation_mode=-1;
end
end
......@@ -359,7 +354,7 @@ if prior_exist % use estimated_params block
if ~isempty(estim_params_.var_exo)
indpstderr = estim_params_.var_exo(:,1); %values correspond to varexo declaration order, row number corresponds to order in estimated_params
end
indpcorr=[]; %initialize matrix for corr paramters
indpcorr=[]; %initialize matrix for corr parameters
if ~isempty(estim_params_.corrx)
indpcorr = estim_params_.corrx(:,1:2); %values correspond to varexo declaration order, row number corresponds to order in estimated_params
end
......@@ -374,11 +369,11 @@ if prior_exist % use estimated_params block
name_tex = cell(totparam_nbr,1); %initialize cell for TeX parameter names
for jj=1:totparam_nbr
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,'$',''); %ordering corresponds to estimated_params
[param_name_temp, param_name_tex_temp]= get_the_name(jj,options_.TeX,M_,estim_params_,options_.varobs);
name_tex{jj,1} =param_name_tex_temp;
name{jj,1} = param_name_temp; %ordering corresponds to estimated_params
else
param_name_temp = get_the_name(jj,options_.TeX,M_,estim_params_,options_);
param_name_temp = get_the_name(jj,options_.TeX,M_,estim_params_,options_.varobs);
name{jj,1} = param_name_temp; %ordering corresponds to estimated_params
end
end
......@@ -392,10 +387,10 @@ else % no estimated_params block, choose all model parameters and all stderr par
totparam_nbr = modparam_nbr+stderrparam_nbr;
name = cellfun(@(x) horzcat('SE_', x), M_.exo_names, 'UniformOutput', false); %names for stderr parameters
name = vertcat(name, M_.param_names);
name_tex = cellfun(@(x) horzcat('$ SE_{', x, '} $'), M_.exo_names, 'UniformOutput', false);
name_tex = vertcat(name_tex, M_.param_names_tex);
name_tex = cellfun(@(x) horzcat('$ SE_{', x, '} $'), M_.exo_names_tex, 'UniformOutput', false);
name_tex = vertcat(name_tex, cellfun(@(x) horzcat('$ ', x, ' $'), M_.param_names_tex, 'UniformOutput', false));
if ~isequal(M_.H,0)
fprintf('\ndynare_identification:: Identification does not support measurement errors (yet) and will ignore them in the following. To test their identifiability, instead define them explicitly as varexo and provide measurement equations in the model definition.\n')
fprintf('\nidentification.run:: Identification does not support measurement errors (yet) and will ignore them in the following. To test their identifiability, instead define them explicitly as varexo and provide measurement equations in the model definition.\n')
end
end
options_ident.name_tex = name_tex;
......@@ -413,13 +408,13 @@ end
% settings dependent on number of parameters
options_ident = set_default_option(options_ident,'max_dim_cova_group',min([2,totparam_nbr-1]));
options_ident.max_dim_cova_group = min([options_ident.max_dim_cova_group,totparam_nbr-1]);
% In brute force search (ident_bruteforce.m) when advanced=1 this option sets the maximum dimension of groups of parameters that best reproduce the behavior of each single model parameter
% In brute force search (identification.bruteforce.m) when advanced=1 this option sets the maximum dimension of groups of parameters that best reproduce the behavior of each single model parameter
options_ident = set_default_option(options_ident,'checks_via_subsets',0);
% 1: uses identification_checks_via_subsets.m to compute problematic parameter combinations
% 0: uses identification_checks.m to compute problematic parameter combinations [default]
% 1: uses identification.checks_via_subsets.m to compute problematic parameter combinations
% 0: uses identification.checks.m to compute problematic parameter combinations [default]
options_ident = set_default_option(options_ident,'max_dim_subsets_groups',min([4,totparam_nbr-1]));
% In identification_checks_via_subsets.m, when checks_via_subsets=1, this option sets the maximum dimension of groups of parameters for which the corresponding rank criteria is checked
% In identification.checks_via_subsets.m, when checks_via_subsets=1, this option sets the maximum dimension of groups of parameters for which the corresponding rank criteria is checked
% store identification options
......@@ -481,8 +476,8 @@ if iload <=0
end
options_ident.tittxt = parameters; %title text for graphs and figures
% perform identification analysis for single point
[ide_moments_point, ide_spectrum_point, ide_minimal_point, ide_hess_point, ide_reducedform_point, ide_dynamic_point, derivatives_info_point, info, error_indicator_point] = ...
identification_analysis(params, indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1); %the 1 at the end implies initialization of persistent variables
[ide_moments_point, ide_spectrum_point, ide_minimal_point, ide_hess_point, ide_reducedform_point, ide_dynamic_point, ~, info, error_indicator_point] = ...
identification.analysis(M_,options_,oo_,bayestopt_,estim_params_,params, indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1); %the 1 at the end implies initialization of persistent variables
if info(1)~=0
% there are errors in the solution algorithm
message = get_error_message(info,options_);
......@@ -495,11 +490,11 @@ if iload <=0
kk=0;
while kk<50 && info(1)
kk=kk+1;
params = Prior.draw();
params = Prior.draw()'; %column vector is expected
options_ident.tittxt = 'Random_prior_params'; %title text for graphs and figures
% perform identification analysis
[ide_moments_point, ide_spectrum_point, ide_minimal_point, ide_hess_point, ide_reducedform_point, ide_dynamic_point, derivatives_info_point, info, error_indicator_point] = ...
identification_analysis(params, indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1);
[ide_moments_point, ide_spectrum_point, ide_minimal_point, ide_hess_point, ide_reducedform_point, ide_dynamic_point, ~, info, error_indicator_point] = ...
identification.analysis(M_,options_,oo_,bayestopt_,estim_params_,params, indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1);
end
end
if info(1)
......@@ -509,11 +504,15 @@ if iload <=0
fprintf('The model did not solve for any of 50 attempts of random samples from the prior\n');
end
fprintf('-----------\n');
warning(orig_warning_state);
return
else
% found a (random) point that solves the model
fprintf('Found a random draw from the priors that solves the model:\n');
disp(params);
labels = name;
headers = {'Name', 'Value'};
lh = cellofchararraymaxlength(labels)+2;
dyntable(options_, 'Feasible draw', headers, labels, params', lh, 10, 6);
fprintf('Identification now continues for this draw.');
parameters = 'Random_prior_params';
parameters_TeX = 'Random prior parameter draw';
......@@ -524,10 +523,11 @@ if iload <=0
save([IdentifDirectoryName '/' fname '_identif.mat'], 'ide_moments_point', 'ide_spectrum_point', 'ide_minimal_point', 'ide_hess_point', 'ide_reducedform_point', 'ide_dynamic_point', 'store_options_ident');
save([IdentifDirectoryName '/' fname '_' parameters '_identif.mat'], 'ide_moments_point', 'ide_spectrum_point', 'ide_minimal_point', 'ide_hess_point', 'ide_reducedform_point', 'ide_dynamic_point', 'store_options_ident');
% display results of identification analysis
disp_identification(params, ide_reducedform_point, ide_moments_point, ide_spectrum_point, ide_minimal_point, name, options_ident);
identification.display(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, parameters_TeX, name_tex);
identification.plot(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
......@@ -539,7 +539,7 @@ if iload <=0
file_index = 0; % initialize counter for files (if options_.MaxNumberOfBytes is reached, we store results in files)
options_MC = options_ident; %store options structure for Monte Carlo analysis
options_MC.advanced = 0; %do not run advanced checking in a Monte Carlo analysis
options_ident.checks_via_subsets = 0; % for Monte Carlo analysis currently only identification_checks and not identification_checks_via_subsets is supported
options_ident.checks_via_subsets = 0; % for Monte Carlo analysis currently only identification.checks and not identification.checks_via_subsets is supported
else
iteration = 1; % iteration equals SampleSize and we are finished
pdraws = []; % to have output object otherwise map_ident.m may crash
......@@ -552,8 +552,8 @@ if iload <=0
end
options_ident.tittxt = []; % clear title text for graphs and figures
% run identification analysis
[ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide_dynamic, ide_derivatives_info, info, error_indicator] = ...
identification_analysis(params, indpmodel, indpstderr, indpcorr, options_MC, dataset_info, prior_exist, 0); % the 0 implies that we do not initialize persistent variables anymore
[ide_moments, ide_spectrum, ide_minimal, ~, ide_reducedform, ide_dynamic, ~, info, error_indicator] = ...
identification.analysis(M_,options_,oo_,bayestopt_,estim_params_,params, indpmodel, indpstderr, indpcorr, options_MC, dataset_info, prior_exist, 0); % the 0 implies that we do not initialize persistent variables anymore
if iteration==0 && info(1)==0 % preallocate storage in the first admissable run
delete([IdentifDirectoryName '/' fname '_identif_*.mat']) % delete previously saved results
......@@ -746,7 +746,7 @@ if iload <=0
end
run_index = 0; % reset index
end
if SampleSize > 1
if SampleSize > 1 && mod(iteration,3)
dyn_waitbar(iteration/SampleSize, h, ['MC identification checks ', int2str(iteration), '/', int2str(SampleSize)]);
end
end
......@@ -796,28 +796,41 @@ if iload <=0
else
maxrun_dMINIMAL = 0;
end
si_dDYNAMICnorm=NaN(max([maxrun_dDYNAMIC, maxrun_dREDUCEDFORM, maxrun_dMOMENTS, maxrun_dSPECTRUM, maxrun_dMINIMAL]),size(STO_si_dDYNAMIC,2));
if ~options_MC.no_identification_reducedform
si_dREDUCEDFORMnorm=NaN(max([maxrun_dDYNAMIC, maxrun_dREDUCEDFORM, maxrun_dMOMENTS, maxrun_dSPECTRUM, maxrun_dMINIMAL]),size(STO_si_dREDUCEDFORM,2));
end
if ~options_MC.no_identification_moments
si_dMOMENTSnorm=NaN(max([maxrun_dDYNAMIC, maxrun_dREDUCEDFORM, maxrun_dMOMENTS, maxrun_dSPECTRUM, maxrun_dMINIMAL]),size(STO_si_dMOMENTS,2));
end
if ~options_MC.no_identification_spectrum
dSPECTRUMnorm=NaN(max([maxrun_dDYNAMIC, maxrun_dREDUCEDFORM, maxrun_dMOMENTS, maxrun_dSPECTRUM, maxrun_dMINIMAL]),size(STO_dSPECTRUM,2));
end
if ~options_MC.no_identification_minimal
dMINIMALnorm=NaN(max([maxrun_dDYNAMIC, maxrun_dREDUCEDFORM, maxrun_dMOMENTS, maxrun_dSPECTRUM, maxrun_dMINIMAL]),size(STO_dMINIMAL,2));
end
for irun=1:max([maxrun_dDYNAMIC, maxrun_dREDUCEDFORM, maxrun_dMOMENTS, maxrun_dSPECTRUM, maxrun_dMINIMAL])
iter=iter+1;
% note that this is not the same si_dDYNAMICnorm as computed in identification_analysis
% note that this is not the same si_dDYNAMICnorm as computed in identification.analysis
% given that we have the MC sample of the Jacobians, we also normalize by the std of the sample of Jacobian entries, to get a fully standardized sensitivity measure
si_dDYNAMICnorm(iter,:) = vnorm(STO_si_dDYNAMIC(:,:,irun)./repmat(normalize_STO_DYNAMIC,1,totparam_nbr-(stderrparam_nbr+corrparam_nbr))).*normaliz1((stderrparam_nbr+corrparam_nbr)+1:end);
si_dDYNAMICnorm(iter,:) = identification.vnorm(STO_si_dDYNAMIC(:,:,irun)./repmat(normalize_STO_DYNAMIC,1,totparam_nbr-(stderrparam_nbr+corrparam_nbr))).*normaliz1((stderrparam_nbr+corrparam_nbr)+1:end);
if ~options_MC.no_identification_reducedform && ~isempty(STO_si_dREDUCEDFORM)
% note that this is not the same si_dREDUCEDFORMnorm as computed in identification_analysis
% note that this is not the same si_dREDUCEDFORMnorm as computed in identification.analysis
% given that we have the MC sample of the Jacobians, we also normalize by the std of the sample of Jacobian entries, to get a fully standardized sensitivity measure
si_dREDUCEDFORMnorm(iter,:) = vnorm(STO_si_dREDUCEDFORM(:,:,irun)./repmat(normalize_STO_REDUCEDFORM,1,totparam_nbr)).*normaliz1;
si_dREDUCEDFORMnorm(iter,:) = identification.vnorm(STO_si_dREDUCEDFORM(:,:,irun)./repmat(normalize_STO_REDUCEDFORM,1,totparam_nbr)).*normaliz1;
end
if ~options_MC.no_identification_moments && ~isempty(STO_si_dMOMENTS)
% note that this is not the same si_dMOMENTSnorm as computed in identification_analysis
% note that this is not the same si_dMOMENTSnorm as computed in identification.analysis
% given that we have the MC sample of the Jacobians, we also normalize by the std of the sample of Jacobian entries, to get a fully standardized sensitivity measure
si_dMOMENTSnorm(iter,:) = vnorm(STO_si_dMOMENTS(:,:,irun)./repmat(normalize_STO_MOMENTS,1,totparam_nbr)).*normaliz1;
si_dMOMENTSnorm(iter,:) = identification.vnorm(STO_si_dMOMENTS(:,:,irun)./repmat(normalize_STO_MOMENTS,1,totparam_nbr)).*normaliz1;
end
if ~options_MC.no_identification_spectrum && ~isempty(STO_dSPECTRUM)
% note that this is not the same dSPECTRUMnorm as computed in identification_analysis
dSPECTRUMnorm(iter,:) = vnorm(STO_dSPECTRUM(:,:,irun)); %not yet used
% note that this is not the same dSPECTRUMnorm as computed in identification.analysis
dSPECTRUMnorm(iter,:) = identification.vnorm(STO_dSPECTRUM(:,:,irun)); %not yet used
end
if ~options_MC.no_identification_minimal && ~isempty(STO_dMINIMAL)
% note that this is not the same dMINIMALnorm as computed in identification_analysis
dMINIMALnorm(iter,:) = vnorm(STO_dMINIMAL(:,:,irun)); %not yet used
% note that this is not the same dMINIMALnorm as computed in identification.analysis
dMINIMALnorm(iter,:) = identification.vnorm(STO_dMINIMAL(:,:,irun)); %not yet used
end
end
end
......@@ -844,7 +857,7 @@ else
options_.options_ident = options_ident;
end
%% if dynare_identification is called as it own function (not through identification command) and if we load files
%% if identification.run is called as it own function (not through identification command) and if we load files
if nargout>3 && iload
filnam = dir([IdentifDirectoryName '/' fname '_identif_*.mat']);
STO_si_dDYNAMIC = [];
......@@ -873,10 +886,11 @@ end
if iload
%if previous analysis is loaded
fprintf(['Testing %s\n',parameters]);
disp_identification(ide_hess_point.params, ide_reducedform_point, ide_moments_point, ide_spectrum_point, ide_minimal_point, name, options_ident);
identification.display(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, [], name_tex);
identification.plot(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
......@@ -886,11 +900,12 @@ if SampleSize > 1
%print results to console but make sure advanced=0
advanced0 = options_ident.advanced;
options_ident.advanced = 0;
disp_identification(pdraws, IDE_REDUCEDFORM, IDE_MOMENTS, IDE_SPECTRUM, IDE_MINIMAL, name, options_ident);
identification.display(pdraws, IDE_REDUCEDFORM, IDE_MOMENTS, IDE_SPECTRUM, IDE_MINIMAL, name, options_ident);
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, [], name_tex);
identification.plot(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
......@@ -906,16 +921,17 @@ if SampleSize > 1
fprintf('\nTesting %s.\n',tittxt);
if ~iload
options_ident.tittxt = tittxt; %title text for graphs and figures
[ide_moments_max, ide_spectrum_max, ide_minimal_max, ide_hess_max, ide_reducedform_max, ide_dynamic_max, derivatives_info_max, info_max, error_indicator_max] = ...
identification_analysis(pdraws(jmax,:), indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1); %the 1 at the end initializes some persistent variables
[ide_moments_max, ide_spectrum_max, ide_minimal_max, ide_hess_max, ide_reducedform_max, ide_dynamic_max, ~, ~, error_indicator_max] = ...
identification.analysis(M_,options_,oo_,bayestopt_,estim_params_,pdraws(jmax,:), indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1); %the 1 at the end initializes some persistent variables
save([IdentifDirectoryName '/' fname '_identif.mat'], 'ide_hess_max', 'ide_moments_max', 'ide_spectrum_max', 'ide_minimal_max','ide_reducedform_max', 'ide_dynamic_max', 'jmax', '-append');
end
advanced0 = options_ident.advanced; options_ident.advanced = 1; % make sure advanced setting is on
disp_identification(pdraws(jmax,:), ide_reducedform_max, ide_moments_max, ide_spectrum_max, ide_minimal_max, name, options_ident);
identification.display(pdraws(jmax,:), ide_reducedform_max, ide_moments_max, ide_spectrum_max, ide_minimal_max, name, options_ident);
options_ident.advanced = advanced0; %reset advanced setting
if ~options_.nograph && ~error_indicator_max.identification_strength && ~error_indicator_max.identification_moments
% plot (i) identification strength and sensitivity measure based on the sample information matrix and (ii) advanced analysis graphs
plot_identification(pdraws(jmax,:), ide_moments_max, ide_hess_max, ide_reducedform_max, ide_dynamic_max, 1, tittxt, name, IdentifDirectoryName, tittxt, name_tex);
identification.plot(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,16 +940,17 @@ if SampleSize > 1
fprintf('Testing %s.\n',tittxt);
if ~iload
options_ident.tittxt = tittxt; %title text for graphs and figures
[ide_moments_min, ide_spectrum_min, ide_minimal_min, ide_hess_min, ide_reducedform_min, ide_dynamic_min, derivatives_info_min, info_min, error_indicator_min] = ...
identification_analysis(pdraws(jmin,:), indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1); %the 1 at the end initializes persistent variables
[ide_moments_min, ide_spectrum_min, ide_minimal_min, ide_hess_min, ide_reducedform_min, ide_dynamic_min, ~, ~, error_indicator_min] = ...
identification.analysis(M_,options_,oo_,bayestopt_,estim_params_,pdraws(jmin,:), indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1); %the 1 at the end initializes persistent variables
save([IdentifDirectoryName '/' fname '_identif.mat'], 'ide_hess_min', 'ide_moments_min','ide_spectrum_min','ide_minimal_min','ide_reducedform_min', 'ide_dynamic_min', 'jmin', '-append');
end
advanced0 = options_ident.advanced; options_ident.advanced = 1; % make sure advanced setting is on
disp_identification(pdraws(jmin,:), ide_reducedform_min, ide_moments_min, ide_spectrum_min, ide_minimal_min, name, options_ident);
identification.display(pdraws(jmin,:), ide_reducedform_min, ide_moments_min, ide_spectrum_min, ide_minimal_min, name, options_ident);
options_ident.advanced = advanced0; %reset advanced setting
if ~options_.nograph && ~error_indicator_min.identification_strength && ~error_indicator_min.identification_moments
% plot (i) identification strength and sensitivity measure based on the sample information matrix and (ii) advanced analysis graphs
plot_identification(pdraws(jmin,:),ide_moments_min,ide_hess_min,ide_reducedform_min,ide_dynamic_min,1,tittxt,name,IdentifDirectoryName,tittxt,name_tex);
identification.plot(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;
......@@ -947,14 +964,15 @@ if SampleSize > 1
if ~iload
options_ident.tittxt = tittxt; %title text for graphs and figures
[ide_moments_(j), ide_spectrum_(j), ide_minimal_(j), ide_hess_(j), ide_reducedform_(j), ide_dynamic_(j), derivatives_info_(j), info_resolve, error_indicator_j] = ...
identification_analysis(pdraws(jcrit(j),:), indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1);
identification.analysis(M_,options_,oo_,bayestopt_,estim_params_,pdraws(jcrit(j),:), indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1);
end
advanced0 = options_ident.advanced; options_ident.advanced = 1; %make sure advanced setting is on
disp_identification(pdraws(jcrit(j),:), ide_reducedform_(j), ide_moments_(j), ide_spectrum_(j), ide_minimal_(j), name, options_ident);
identification.display(pdraws(jcrit(j),:), ide_reducedform_(j), ide_moments_(j), ide_spectrum_(j), ide_minimal_(j), name, options_ident);
options_ident.advanced = advanced0; % reset advanced
if ~options_.nograph && ~error_indicator_j.identification_strength && ~error_indicator_j.identification_moments
% plot (i) identification strength and sensitivity measure based on the sample information matrix and (ii) advanced analysis graphs
plot_identification(pdraws(jcrit(j),:), ide_moments_(j), ide_hess_(j), ide_reducedform_(j), ide_dynamic_(j), 1, tittxt, name, IdentifDirectoryName, tittxt, name_tex);
identification.plot(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
......@@ -968,8 +986,6 @@ if SampleSize > 1
end
%reset warning state
warning_config;
warning(orig_warning_state);
fprintf('\n==== Identification analysis completed ====\n\n')
options_ = store_options_; %restore options set
matlab/simulated_moment_uncertainty.m → matlab/+identification/simulated_moment_uncertainty.m
View file @ 4a32acc2
......@@ -12,7 +12,7 @@ function [cmm, mm] = simulated_moment_uncertainty(indx, periods, replic,options_
% - cmm: [n_moments by n_moments] covariance matrix of simulated moments
% - mm: [n_moments by replic] matrix of moments
% Copyright © 2009-2018 Dynare Team
% Copyright © 2009-2024 Dynare Team
%
% This file is part of Dynare.
%
......@@ -58,17 +58,12 @@ if M_.exo_nbr > 0
oo_.exo_simul= ones(max(options_.periods,1) + M_.maximum_lag + M_.maximum_lead,1) * oo_.exo_steady_state';
end
oo_.dr=set_state_space(oo_.dr,M_,options_);
oo_.dr=set_state_space(oo_.dr,M_);
if options_.logged_steady_state %if steady state was previously logged, undo this
oo_.dr.ys=exp(oo_.dr.ys);
oo_.steady_state=exp(oo_.steady_state);
options_.logged_steady_state=0;
logged_steady_state_indicator=1;
evalin('base','options_.logged_steady_state=0;')
else
logged_steady_state_indicator=0;
end
[oo_.dr,info,M_.params] = compute_decision_rules(M_,options_,oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
......@@ -92,7 +87,6 @@ else
end
end
for j=1:replic
[ys, oo_.exo_simul] = simult(y0,oo_.dr,M_,options_);%do simulation
oo_=disp_moments(ys, options_.varobs,M_,options_,oo_); %get moments
......@@ -106,8 +100,5 @@ for j=1:replic
end
dyn_waitbar_close(h);
if logged_steady_state_indicator
evalin('base','options_.logged_steady_state=1;') %reset base workspace option to conform to base oo_
end
cmm = cov(mm');
disp('Simulated moment uncertainty ... done!')
matlab/unfold_g3.m → matlab/+identification/unfold_g3.m
View file @ 4a32acc2
function g3_unfolded = unfold_g3(g3, ny)
% Given the 3rd order derivatives stored in a sparse matrix and without
% symmetric elements (as returned by the static/dynamic files) and the number
% of (static or dynamic )variables in the jacobian, returns
% of (static or dynamic )variables in the Jacobian, returns
% an unfolded version of the same matrix (i.e. with symmetric elements).
% Copyright © 2019 Dynare Team
......
matlab/unfold_g4.m → matlab/+identification/unfold_g4.m
View file @ 4a32acc2
function g4_unfolded = unfold_g4(g4, ny)
% Given the 4th order derivatives stored in a sparse matrix and without
% symmetric elements (as returned by the static/dynamic files) and the number
% of (static or dynamic) variables in the jacobian, returns
% of (static or dynamic) variables in the Jacobian, returns
% an unfolded version of the same matrix (i.e. with symmetric elements).
% Copyright © 2019 Dynare Team
......
matlab/vnorm.m → matlab/+identification/vnorm.m
View file @ 4a32acc2
File moved
matlab/+mom/Jtest.m
View file @ 4a32acc2
function oo_ = Jtest(xparam, objective_function, Woptflag, oo_, options_mom_, bayestopt_, Bounds, estim_params_, M_, nobs)
% function oo_ = Jtest(xparam, objective_function, Woptflag, oo_, options_mom_, bayestopt_, Bounds, estim_params_, M_, nobs)
function J_test = Jtest(xparam, objective_function, Q, model_moments, m_data, data_moments, weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
% J_test = Jtest(xparam, objective_function, Q, model_moments, m_data, data_moments, weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
% -------------------------------------------------------------------------
% Computes the J-test statistic and p-value given a GMM/SMM estimation
% -------------------------------------------------------------------------
% Computes the J-test statistic and p-value for a GMM/SMM estimation
% =========================================================================
% INPUTS
% xparam: [vector] estimated parameter vector
% objective_function: [function handle] objective function
% Woptflag: [logical] flag if optimal weighting matrix has already been computed
% oo_: [struct] results
% Q: [scalar] value of moments distance criterion
% model_moments: [vector] model moments
% m_data: [matrix] selected empirical moments at each point in time
% data_moments: [vector] empirical moments
% weighting_info: [struct] information on weighting matrix
% options_mom_: [struct] options
% bayestopt_: [struct] information on priors
% Bounds: [struct] bounds on parameters
% estim_params_: [struct] information on estimated parameters
% M_: [struct] information on the model
% nobs: [scalar] number of observations
% M_: [struct] model information
% estim_params_: [struct] estimated parameters
% bayestopt_: [struct] info on prior distributions
% BoundsInfo: [struct] info bounds on parameters
% dr: [struct] reduced form model
% endo_steady_state: [vector] steady state of endogenous variables (initval)
% exo_steady_state: [vector] steady state of exogenous variables (initval)
% exo_det_steady_state: [vector] steady state of deterministic exogenous variables (initval)
% -------------------------------------------------------------------------
% OUTPUT
% oo_: [struct] updated results
% J_test: [struct] results of J test
% -------------------------------------------------------------------------
% This function is called by
% o mom.run
......@@ -24,7 +30,8 @@ function oo_ = Jtest(xparam, objective_function, Woptflag, oo_, options_mom_, ba
% This function calls
% o mom.objective_function
% o mom.optimal_weighting_matrix
% =========================================================================
% -------------------------------------------------------------------------
% Copyright © 2023 Dynare Team
%
% This file is part of Dynare.
......@@ -41,27 +48,28 @@ function oo_ = Jtest(xparam, objective_function, Woptflag, oo_, options_mom_, ba
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% =========================================================================
if options_mom_.mom.mom_nbr > length(xparam)
% Get optimal weighting matrix for J test, if necessary
if ~Woptflag
W_opt = mom.optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.model_moments, options_mom_.mom.bartlett_kernel_lag);
oo_J = oo_;
oo_J.mom.Sw = chol(W_opt);
fval = feval(objective_function, xparam, Bounds, oo_J, estim_params_, M_, options_mom_, bayestopt_);
if ~weighting_info.Woptflag
W_opt = mom.optimal_weighting_matrix(m_data, model_moments, options_mom_.mom.bartlett_kernel_lag);
weighting_info.Sw = chol(W_opt);
fval = feval(objective_function, xparam, data_moments, weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
else
fval = oo_.mom.Q;
fval = Q;
end
% Compute J statistic
if strcmp(options_mom_.mom.mom_method,'SMM')
Variance_correction_factor = options_mom_.mom.variance_correction_factor;
variance_correction_factor = options_mom_.mom.variance_correction_factor;
elseif strcmp(options_mom_.mom.mom_method,'GMM')
Variance_correction_factor = 1;
variance_correction_factor = 1;
end
oo_.mom.J_test.j_stat = nobs*Variance_correction_factor*fval/options_mom_.mom.weighting_matrix_scaling_factor;
oo_.mom.J_test.degrees_freedom = length(oo_.mom.model_moments)-length(xparam);
oo_.mom.J_test.p_val = 1-chi2cdf(oo_.mom.J_test.j_stat, oo_.mom.J_test.degrees_freedom);
fprintf('\nValue of J-test statistic: %f\n',oo_.mom.J_test.j_stat);
fprintf('p-value of J-test statistic: %f\n',oo_.mom.J_test.p_val);
J_test.j_stat = options_mom_.nobs*variance_correction_factor*fval/options_mom_.mom.weighting_matrix_scaling_factor;
J_test.degrees_freedom = length(model_moments)-length(xparam);
J_test.p_val = 1-chi2cdf(J_test.j_stat, J_test.degrees_freedom);
fprintf('\nValue of J-test statistic: %f\n',J_test.j_stat);
fprintf('p-value of J-test statistic: %f\n',J_test.p_val);
else
J_test=[];
end
\ No newline at end of file
matlab/+mom/check_irf_matching_file.m 0 → 100644
View file @ 4a32acc2
function [irf_matching_file_name, irf_matching_file_path] = check_irf_matching_file(irf_matching_file)
% [irf_matching_file_name, irf_matching_file_path] = check_irf_matching_file(irf_matching_file)
% -------------------------------------------------------------------------
% Check if the provided irf_matching_file is a valid MATLAB function with
% .m extension and return name, path and extension of the file.
% -------------------------------------------------------------------------
% INPUTS
% - irf_matching_file: [string] user provided name (with possible path and extension)
% of the MATLAB function that transforms model IRFs
% -------------------------------------------------------------------------
% OUTPUTS
% - irf_matching_file_name: [string] name of the MATLAB function (without extension)
% - irf_matching_file_path: [string] path to irf_matching_file_name
% -------------------------------------------------------------------------
% This function is called by
% - mom.run
% -------------------------------------------------------------------------
% Copyright © 2023 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
if isempty(irf_matching_file)
% no irf_matching_file provided, so no transformations will be done
irf_matching_file_name = '';
irf_matching_file_path = '';
else
[irf_matching_file_path, irf_matching_file_name, irf_matching_file_ext] = fileparts(irf_matching_file);
% make sure file is a MATLAB function with .m extension
if ~strcmp(irf_matching_file_ext,'.m')
if strcmp(irf_matching_file_ext,'')
irf_matching_file_ext = '.m';
else
error('method_of_moments: ''irf_matching_file'' needs to point towards a MATLAB function with extension ''.m''!');
end
end
if isempty(irf_matching_file_path)
irf_matching_file_path = '.';
end
if exist([irf_matching_file_path filesep irf_matching_file_name irf_matching_file_ext],'file') ~= 2
error('method_of_moments: Could not find a ''irf_matching_file'' called ''%s''!',[irf_matching_file_path filesep irf_matching_file_name irf_matching_file_ext]);
end
end
\ No newline at end of file
matlab/+mom/default_option_mom_values.m
View file @ 4a32acc2
function options_mom_ = default_option_mom_values(options_mom_, options_, dname, doBayesianEstimation)
% function options_mom_ = default_option_mom_values(options_mom_, options_, dname, doBayesianEstimation)
% Returns structure containing the options for method_of_moments command
% options_mom_ is local and contains default and user-specified values for
% all settings needed for the method of moments estimation. Some options,
% though, are set by the preprocessor into options_ and we copy these over.
% The idea is to be independent of options_ and have full control of the
% estimation instead of possibly having to deal with options chosen somewhere
% else in the mod file.
% =========================================================================
function options_mom_ = default_option_mom_values(options_mom_, options_, dname, fname, do_bayesian_estimation)
% options_mom_ = default_option_mom_values(options_mom_, options_, dname, fname, do_bayesian_estimation)
% -------------------------------------------------------------------------
% Returns structure containing the options for method_of_moments command.
% Note 1: options_mom_ is local and contains default and user-specified
% values for all settings needed for the method of moments estimation.
% Some options, though, are set by the preprocessor into options_ and we
% copy these over. The idea is to be independent of options_ and have full
% control of the estimation instead of possibly having to deal with options
% chosen somewhere else in the mod file.
% Note 2: we call a "mode" the minimum of the objective function, i.e.
% the parameter vector that minimizes the distance between the moments/IRFs
% computed from the model and the moments/IRFs computed from the data.
% -------------------------------------------------------------------------
% INPUTS
% o options_mom_: [structure] information about all (user-specified and updated) settings used in estimation (options_mom_)
% o options_: [structure] information on global options
% o dname: [string] name of directory to store results
% o doBayesianEstimation [boolean] indicator whether we do Bayesian estimation
% o options_mom_: [structure] all user-specified settings (from the method_of_moments command)
% o options_: [structure] global options
% o dname: [string] default name of directory to store results
% o fname: [string] default name of mod file
% o do_bayesian_estimation [boolean] indicator whether we do Bayesian estimation
% -------------------------------------------------------------------------
% OUTPUTS
% o oo_: [structure] storage for results (oo_)
% o options_mom_: [structure] information about all (user-specified and updated) settings used in estimation (options_mom_)
% o options_mom_: [structure] all user-specified and updated settings required for method_of_moments estimation
% -------------------------------------------------------------------------
% This function is called by
% o mom.run
......@@ -46,34 +47,43 @@ function options_mom_ = default_option_mom_values(options_mom_, options_, dname,
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% =========================================================================
mom_method = options_mom_.mom.mom_method; % this is a required option
% -------------------------------------------------------------------------
% LIMITATIONS
% -------------------------------------------------------------------------
if options_.logged_steady_state || options_.loglinear
error('method_of_moments: The loglinear option is not supported. Please append the required logged variables as auxiliary equations.')
error('method_of_moments: The loglinear option is not supported. Please append the required logged variables as auxiliary equations.');
else
options_mom_.logged_steady_state = 0;
options_mom_.loglinear = false;
end
options_mom_.hessian.use_penalized_objective = false; % penalized objective not yet
% options related to variable declarations
if isfield(options_mom_,'hessian') && options_mom_.hessian.use_penalized_objective
warning('method_of_moments: The ''use_penalized_objective_for_hessian'' option is not supported yet and will be skipped.');
end
options_mom_.hessian.use_penalized_objective = false; % penalized objective not yet supported
if isfield(options_,'trend_coeffs')
error('method_of_moments: %s does not allow for trend in data',mom_method)
error('method_of_moments: %s does not allow for trend in data',mom_method);
end
% options related to endogenous prior restrictions are not supported
if ~isempty(options_.endogenous_prior_restrictions.irf) && ~isempty(options_.endogenous_prior_restrictions.moment)
fprintf('method_of_moments: Endogenous prior restrictions are not supported yet and will be skipped.\n')
warning('method_of_moments: Endogenous prior restrictions are not supported yet and will be skipped.');
end
options_mom_.endogenous_prior_restrictions.irf = {};
options_mom_.endogenous_prior_restrictions.moment = {};
options_mom_.mom.analytic_jacobian_optimizers = [1, 3, 4, 13, 101]; % these are currently supported optimizers that are able to use the analytical_jacobian option
if isfield(options_mom_,'bayesian_irf') && options_mom_.bayesian_irf % do we need this at all??
warning('method_of_moments: The ''bayesian_irf'' option is not supported yet and will be skipped.');
end
options_mom_.bayesian_irf = false;
if strcmp(mom_method,'IRF_MATCHING')
if isfield(options_mom_.mom,'penalized_estimator') && options_mom_.mom.penalized_estimator
warning('method_of_moments: The ''penalized_estimator'' option is not supported yet for IRF_MATCHING and will be ignored.');
end
options_mom_.mom.penalized_estimator = false;
end
% -------------------------------------------------------------------------
% OPTIONS POSSIBLY SET BY THE USER
......@@ -87,8 +97,7 @@ options_mom_ = set_default_option(options_mom_,'nograph',false); % do no
options_mom_ = set_default_option(options_mom_,'noprint',false); % do not print output to console
options_mom_ = set_default_option(options_mom_,'TeX',false); % print TeX tables and graphics
options_mom_.mom = set_default_option(options_mom_.mom,'verbose',false); % display and store intermediate estimation results
%options_mom_ = set_default_option(options_mom_,'verbosity',false); %
if doBayesianEstimation
if do_bayesian_estimation
options_mom_ = set_default_option(options_mom_,'plot_priors',true); % control plotting of priors
options_mom_ = set_default_option(options_mom_,'prior_trunc',1e-10); % probability of extreme values of the prior density that is ignored when computing bounds for the parameters
end
......@@ -111,6 +120,16 @@ end
if strcmp(mom_method,'GMM')
options_mom_.mom = set_default_option(options_mom_.mom,'analytic_standard_errors',false); % compute standard errors numerically (0) or analytically (1). Analytical derivatives are only available for GMM.
end
if strcmp(mom_method,'IRF_MATCHING')
if ~isfield(options_mom_.mom,'irf_matching_file')
options_mom_.mom.irf_matching_file = []; % irf_matching file enables to transform model IRFs before matching them to data IRFs
end
options_mom_.mom.irf_matching_file = set_default_option(options_mom_.mom.irf_matching_file,'name','');
options_mom_.mom = set_default_option(options_mom_.mom,'simulation_method','STOCH_SIMUL'); % simulation method used to compute IRFs
options_mom_ = set_default_option(options_mom_,'add_tiny_number_to_cholesky',1e-14); % add tiny number to Cholesky factor to avoid numerical problems when computing IRFs
options_mom_ = set_default_option(options_mom_,'drop',100); % truncation / burnin for order>1 irf simulations
options_mom_ = set_default_option(options_mom_,'relative_irf',false); % requests the computation of normalized IRFs
end
% data related options
if strcmp(mom_method,'GMM') || strcmp(mom_method,'SMM')
......@@ -123,12 +142,14 @@ if strcmp(mom_method,'GMM') || strcmp(mom_method,'SMM')
end
% optimization related
if (isoctave && user_has_octave_forge_package('optim')) || (~isoctave && user_has_matlab_license('optimization_toolbox'))
if strcmp(mom_method,'GMM') || strcmp(mom_method,'SMM')
if (isoctave && user_has_octave_forge_package('optim')) || (~isoctave && user_has_matlab_license('optimization_toolbox'))
options_mom_ = set_default_option(options_mom_,'mode_compute',13); % specifies lsqnonlin as default optimizer for minimization
end
else
options_mom_ = set_default_option(options_mom_,'mode_compute',4); % specifies csminwel as fallback default option for minimization
options_mom_ = set_default_option(options_mom_,'mode_compute',5); % specifies newrat as fallback default option for minimization
end
elseif strcmp(mom_method,'IRF_MATCHING')
options_mom_ = set_default_option(options_mom_,'mode_compute',5); % specifies newrat as fallback default option for minimization
end
options_mom_ = set_default_option(options_mom_,'additional_optimizer_steps',[]); % vector of additional mode-finders run after mode_compute
options_mom_ = set_default_option(options_mom_,'optim_opt',[]); % a list of NAME and VALUE pairs to set options for the optimization routines. Available options depend on mode_compute
......@@ -136,17 +157,29 @@ options_mom_ = set_default_option(options_mom_,'silent_optimizer',false);
options_mom_ = set_default_option(options_mom_,'huge_number',1e7); % value for replacing the infinite bounds on parameters by finite numbers. Used by some optimizers for numerical reasons
options_mom_.mom = set_default_option(options_mom_.mom,'analytic_jacobian',false); % use analytic Jacobian in optimization, only available for GMM and gradient-based optimizers
options_mom_.optimizer_vec = [options_mom_.mode_compute;num2cell(options_mom_.additional_optimizer_steps)];
options_mom_.mom.analytic_jacobian_optimizers = [1, 3, 4, 13, 101]; % these are currently supported optimizers that are able to use the analytic_jacobian option
options_mom_.analytic_derivation = 0; % force to 0 as we check this seperately in dynare_minimize_objective.m
options_mom_ = set_default_option(options_mom_,'mode_file',''); % name of the file containing initial values for the mode
options_mom_ = set_default_option(options_mom_,'cova_compute',true); % 1: computed covariance via Hessian after the computation of the mode, 0: turn off computation of covariance matrix
% perturbation related
options_mom_ = set_default_option(options_mom_,'order',1); % order of Taylor approximation in perturbation
if strcmp(mom_method,'IRF_MATCHING') % number of simulated series used to compute IRFs
if options_mom_.order == 1
options_mom_ = set_default_option(options_mom_,'replic',1);
else
options_mom_ = set_default_option(options_mom_,'replic',50);
end
end
options_mom_ = set_default_option(options_mom_,'pruning',false); % use pruned state space system at order>1
options_mom_ = set_default_option(options_mom_,'aim_solver',false); % use AIM algorithm to compute perturbation approximation instead of mjdgges
options_mom_ = set_default_option(options_mom_,'k_order_solver',false); % use k_order_perturbation instead of mjdgges
options_mom_ = set_default_option(options_mom_,'dr_cycle_reduction',false); % use cycle reduction algorithm to solve the polynomial equation for retrieving the coefficients associated to the endogenous variables in the decision rule
options_mom_ = set_default_option(options_mom_,'dr_cycle_reduction_maxiter',100); % maximum number of iterations used in the cycle reduction algorithm
options_mom_ = set_default_option(options_mom_,'dr_cycle_reduction_tol',1e-7); % convergence criterion used in the cycle reduction algorithm
options_mom_ = set_default_option(options_mom_,'dr_logarithmic_reduction',false); % use logarithmic reduction algorithm to solve the polynomial equation for retrieving the coefficients associated to the endogenous variables in the decision rule
options_mom_ = set_default_option(options_mom_,'dr_logarithmic_reduction_maxiter',100); % maximum number of iterations used in the logarithmic reduction algorithm
options_mom_ = set_default_option(options_mom_,'dr_logarithmic_reduction_tol',1e-12); % convergence criterion used in the cycle reduction algorithm
options_mom_ = set_default_option(options_mom_,'dr_logarithmic_reduction_tol',1e-12); % convergence criterion used in the logarithmic reduction algorithm
options_mom_ = set_default_option(options_mom_,'qz_criterium',1-1e-6); % value used to split stable from unstable eigenvalues in reordering the Generalized Schur decomposition used for solving first order problems
% if there are no unit roots one can use 1.0 (or slightly below) which we set as default; if they are possible, you may have have multiple unit roots and the accuracy decreases when computing the eigenvalues in lyapunov_symm
% Note that unit roots are only possible at first-order, at higher order we set it to 1 in pruned_state_space_system and focus only on stationary observables.
......@@ -160,15 +193,122 @@ options_mom_ = set_default_option(options_mom_,'lyapunov_srs',false);
options_mom_ = set_default_option(options_mom_,'lyapunov_complex_threshold',1e-15); % complex block threshold for the upper triangular matrix in symmetric Lyapunov equation solver
options_mom_ = set_default_option(options_mom_,'lyapunov_fixed_point_tol',1e-10); % convergence criterion used in the fixed point Lyapunov solver
options_mom_ = set_default_option(options_mom_,'lyapunov_doubling_tol',1e-16); % convergence criterion used in the doubling algorithm
options_mom_ = set_default_option(options_mom_,'sylvester_fp',false); % determines whether to use fixed point algorihtm to solve Sylvester equation (gensylv_fp), faster for large scale models
options_mom_ = set_default_option(options_mom_,'sylvester_fixed_point_tol',1e-12); % convergence criterion used in the fixed point Sylvester solver
% mode check plot
% Bayesian MCMC related
if do_bayesian_estimation
options_mom_ = set_default_option(options_mom_,'mh_replic',0); % number of draws in Metropolis-Hastings and slice samplers
options_mom_ = set_default_option(options_mom_,'mh_posterior_mode_estimation',false); % skip optimizer-based mode-finding and instead compute the mode based on a run of a MCMC
options_mom_ = set_default_option(options_mom_,'load_mh_file',false); % add to previous Metropolis-Hastings or slice simulations instead of starting from scratch
options_mom_ = set_default_option(options_mom_,'load_results_after_load_mh',false); % load the previously computed convergence diagnostics, marginal data density, and posterior statistics from an existing mom_results file instead of recomputing them
if options_mom_.mh_replic > 0 || options_mom_.load_mh_file
options_mom_ = set_default_option(options_mom_,'sub_draws',[]);
options_mom_ = set_default_option(options_mom_,'posterior_max_subsample_draws',1200);
options_mom_ = set_default_option(options_mom_,'mh_nblck',2); % number of parallel chains for Metropolis-Hastings or slice algorithm
options_mom_ = set_default_option(options_mom_,'mh_drop',0.5); % fraction of initially generated parameter vectors to be dropped as a burn-in before using posterior simulations
options_mom_ = set_default_option(options_mom_,'mh_conf_sig',0.9); % confidence/HPD interval used for the computation of prior and posterior statistics
options_mom_ = set_default_option(options_mom_,'mh_recover',false); % attempts to recover a Metropolis-Hastings simulation that crashed prematurely
options_mom_ = set_default_option(options_mom_,'MCMC_jumping_covariance','hessian'); % which covariance to use for the proposal density of the MCMC sampler
if ~isfield(options_mom_,'mh_initialize_from_previous_mcmc')
options_mom_.mh_initialize_from_previous_mcmc.status = false; % pick initial values for new MCMC from a previous one
end
options_mom_.mh_initialize_from_previous_mcmc = set_default_option(options_mom_.mh_initialize_from_previous_mcmc,'directory',''); % pick initial values for new MCMC from a previous one: directory
options_mom_.mh_initialize_from_previous_mcmc = set_default_option(options_mom_.mh_initialize_from_previous_mcmc,'record',''); % pick initial values for new MCMC from a previous one: record file name
options_mom_.mh_initialize_from_previous_mcmc = set_default_option(options_mom_.mh_initialize_from_previous_mcmc,'prior',''); % pick initial values for new MCMC from a previous one: prior file name
if ~isfield(options_mom_,'posterior_sampler_options')
options_mom_.posterior_sampler_options = [];
end
options_mom_.posterior_sampler_options = set_default_option(options_mom_.posterior_sampler_options,'posterior_sampling_method','random_walk_metropolis_hastings'); % selects the sampler used to sample from the posterior distribution during Bayesian estimation
options_mom_.posterior_sampler_options = set_default_option(options_mom_.posterior_sampler_options,'sampling_opt',[]); % used to set options for the posterior sampling methods
switch options_mom_.posterior_sampler_options.posterior_sampling_method
case 'random_walk_metropolis_hastings'
if ~isfield(options_mom_.posterior_sampler_options,'rwmh')
options_mom_.posterior_sampler_options.rwmh = [];
end
options_mom_.posterior_sampler_options.rwmh = set_default_option(options_mom_.posterior_sampler_options.rwmh,'proposal_distribution','rand_multivariate_normal');
options_mom_.posterior_sampler_options.rwmh = set_default_option(options_mom_.posterior_sampler_options.rwmh,'student_degrees_of_freedom',3);
options_mom_.posterior_sampler_options.rwmh = set_default_option(options_mom_.posterior_sampler_options.rwmh,'use_mh_covariance_matrix',false);
options_mom_.posterior_sampler_options.rwmh = set_default_option(options_mom_.posterior_sampler_options.rwmh,'save_tmp_file',false);
case 'tailored_random_block_metropolis_hastings'
if ~isfield(options_mom_.posterior_sampler_options,'tarb')
options_mom_.posterior_sampler_options.tarb = [];
end
options_mom_.posterior_sampler_options.tarb = set_default_option(options_mom_.posterior_sampler_options.tarb,'proposal_distribution','rand_multivariate_normal');
options_mom_.posterior_sampler_options.tarb = set_default_option(options_mom_.posterior_sampler_options.tarb,'student_degrees_of_freedom',3);
options_mom_.posterior_sampler_options.tarb = set_default_option(options_mom_.posterior_sampler_options.tarb,'mode_compute',4);
options_mom_.posterior_sampler_options.tarb = set_default_option(options_mom_.posterior_sampler_options.tarb,'new_block_probability',0.25);
options_mom_.posterior_sampler_options.tarb = set_default_option(options_mom_.posterior_sampler_options.tarb,'optim_opt','');
options_mom_.posterior_sampler_options.tarb = set_default_option(options_mom_.posterior_sampler_options.tarb,'save_tmp_file',true);
case 'slice'
if ~isfield(options_mom_.posterior_sampler_options,'slice')
options_mom_.posterior_sampler_options.slice = [];
end
options_mom_.posterior_sampler_options.slice = set_default_option(options_mom_.posterior_sampler_options.slice,'proposal_distribution','');
options_mom_.posterior_sampler_options.slice = set_default_option(options_mom_.posterior_sampler_options.slice,'rotated',0);
options_mom_.posterior_sampler_options.slice = set_default_option(options_mom_.posterior_sampler_options.slice,'slice_initialize_with_mode',false); % must be used with rotated
options_mom_.posterior_sampler_options.slice = set_default_option(options_mom_.posterior_sampler_options.slice,'use_mh_covariance_matrix',false); % must be used with rotated
options_mom_.posterior_sampler_options.slice = set_default_option(options_mom_.posterior_sampler_options.slice,'WR',[]);
options_mom_.posterior_sampler_options.slice = set_default_option(options_mom_.posterior_sampler_options.slice,'mode_files',[]);
options_mom_.posterior_sampler_options.slice = set_default_option(options_mom_.posterior_sampler_options.slice,'mode',[]);
options_mom_.posterior_sampler_options.slice = set_default_option(options_mom_.posterior_sampler_options.slice,'initial_step_size',0.8);
options_mom_.posterior_sampler_options.slice = set_default_option(options_mom_.posterior_sampler_options.slice,'save_tmp_file',true);
case 'independent_metropolis_hastings'
if ~isfield(options_mom_.posterior_sampler_options,'imh')
options_mom_.posterior_sampler_options.imh = [];
end
options_mom_.posterior_sampler_options.imh = set_default_option(options_mom_.posterior_sampler_options.imh,'proposal_distribution','rand_multivariate_normal');
options_mom_.posterior_sampler_options.imh = set_default_option(options_mom_.posterior_sampler_options.imh,'use_mh_covariance_matrix',false);
options_mom_.posterior_sampler_options.imh = set_default_option(options_mom_.posterior_sampler_options.imh,'save_tmp_file',false);
end
if ~strcmp(options_mom_.posterior_sampler_options.posterior_sampling_method,'slice')
options_mom_ = set_default_option(options_mom_,'mh_init_scale_factor',2);
options_mom_ = set_default_option(options_mom_,'mh_jscale',[]);
end
% mh_tune_jscale options
if strcmp(options_mom_.posterior_sampler_options.posterior_sampling_method,'random_walk_metropolis_hastings')
if ~isfield(options_mom_,'mh_tune_jscale')
options_mom_.mh_tune_jscale = [];
end
options_mom_.mh_tune_jscale = set_default_option(options_mom_.mh_tune_jscale,'status',false);
options_mom_.mh_tune_jscale = set_default_option(options_mom_.mh_tune_jscale,'target',0.33);
options_mom_.mh_tune_jscale = set_default_option(options_mom_.mh_tune_jscale,'guess',[]);
options_mom_.mh_tune_jscale.maxiter = options_.mh_tune_jscale.maxiter;
options_mom_.mh_tune_jscale.rho = options_.mh_tune_jscale.rho;
options_mom_.mh_tune_jscale.stepsize = options_.mh_tune_jscale.stepsize;
options_mom_.mh_tune_jscale.c1 = options_.mh_tune_jscale.c1;
options_mom_.mh_tune_jscale.c2 = options_.mh_tune_jscale.c2;
options_mom_.mh_tune_jscale.c3 = options_.mh_tune_jscale.c3;
end
% convergence diagnostics
options_mom_ = set_default_option(options_mom_,'nodiagnostic',false);
if ~isfield(options_mom_,'convergence')
options_mom_.convergence = [];
end
if ~isfield(options_mom_.convergence,'geweke')
options_mom_.convergence.geweke = [];
end
if ~isfield(options_mom_.convergence,'rafterylewis')
options_mom_.convergence.rafterylewis = [];
end
if ~isfield(options_mom_.convergence,'brooksgelman')
options_mom_.convergence.brooksgelman = [];
end
options_mom_.convergence.geweke = set_default_option(options_mom_.convergence.geweke,'taper_steps', [4 8 15]);
options_mom_.convergence.geweke = set_default_option(options_mom_.convergence.geweke,'geweke_interval', [0.2 0.5]);
options_mom_.convergence.rafterylewis = set_default_option(options_mom_.convergence.rafterylewis,'indicator', false);
options_mom_.convergence.rafterylewis = set_default_option(options_mom_.convergence.rafterylewis,'qrs', [0.025 0.005 0.95]);
options_mom_.convergence.brooksgelman = set_default_option(options_mom_.convergence.brooksgelman,'plotrows',3);
end
end
% mode check plot options
options_mom_.mode_check.nolik = false; % we don't do likelihood (also this initializes mode_check substructure)
options_mom_.mode_check = set_default_option(options_mom_.mode_check,'status',false); % plot the target function for values around the computed minimum for each estimated parameter in turn. This is helpful to diagnose problems with the optimizer.
options_mom_.mode_check = set_default_option(options_mom_.mode_check,'neighbourhood_size',.5); % width of the window around the computed minimum to be displayed on the diagnostic plots. This width is expressed in percentage deviation. The Inf value is allowed, and will trigger a plot over the entire domain
options_mom_.mode_check = set_default_option(options_mom_.mode_check,'symmetric_plots',true); % ensure that the check plots are symmetric around the minimum. A value of 0 allows to have asymmetric plots, which can be useful if the minimum is close to a domain boundary, or in conjunction with neighbourhood_size = Inf when the domain is not the entire real line
options_mom_.mode_check = set_default_option(options_mom_.mode_check,'number_of_points',20); % number of points around the minimum where the target function is evaluated (for each parameter)
options_mom_.mode_check = set_default_option(options_mom_.mode_check,'status',false); % plot the target function for values around the computed mode for each estimated parameter in turn. This is helpful to diagnose problems with the optimizer.
options_mom_.mode_check = set_default_option(options_mom_.mode_check,'neighbourhood_size',.5); % width of the window around the computed mode to be displayed on the diagnostic plots. This width is expressed in percentage deviation. The Inf value is allowed, and will trigger a plot over the entire domain
options_mom_.mode_check = set_default_option(options_mom_.mode_check,'symmetric_plots',true); % ensure that the check plots are symmetric around the mode. A value of 0 allows to have asymmetric plots, which can be useful if the mode is close to a domain boundary, or in conjunction with neighbourhood_size = Inf when the domain is not the entire real line
options_mom_.mode_check = set_default_option(options_mom_.mode_check,'number_of_points',20); % number of points around the mode where the target function is evaluated (for each parameter)
% -------------------------------------------------------------------------
......@@ -180,10 +320,19 @@ options_mom_.varobs = options_.varobs; % observable variables in or
options_mom_.varobs_id = options_.varobs_id; % index for observable variables in M_.endo_names
options_mom_.obs_nbr = length(options_mom_.varobs); % number of observed variables
%related to discretionary policy
options_mom_.discretionary_policy = options_.discretionary_policy; % whether discretionary policy is used
options_mom_.discretionary_tol = options_.discretionary_tol; % tolerance of discretionary policy
options_mom_.dp = options_.dp; % other options
% related to call of dynare
options_mom_.console_mode = options_.console_mode;
if options_mom_.console_mode
options_mom_.nodisplay = true;
end
options_mom_.parallel = options_.parallel;
options_mom_.parallel_info = options_.parallel_info;
options_mom_.debug = options_.debug; % debug option is needed by some functions, e.g. check_plot
% related to estimated_params and estimated_params_init blocks
options_mom_.use_calibration_initialization = options_.use_calibration_initialization;
......@@ -252,19 +401,82 @@ end
options_mom_.gstep = options_.gstep; % needed by hessian.m
options_mom_.trust_region_initial_step_bound_factor = options_.trust_region_initial_step_bound_factor; % used in dynare_solve for trust_region
% other
% miscellaneous
options_mom_.threads = options_.threads;
options_mom_.MaxNumberOfBytes = options_.MaxNumberOfBytes;
%options_mom_.MaximumNumberOfMegaBytes = options_.MaximumNumberOfMegaBytes;
options_mom_.marginal_data_density = options_.marginal_data_density;
% -------------------------------------------------------------------------
% DEFAULT VALUES
% -------------------------------------------------------------------------
options_mom_.analytic_derivation = 0;
options_mom_.mom.compute_derivs = false; % flag to compute derivs in objective function (might change for GMM with either analytic_standard_errors or analytic_jacobian (dependent on optimizer))
options_mom_.mom.vector_output = false; % specifies whether the objective function returns a vector
options_mom_.analytic_derivation_mode = 0; % needed by get_perturbation_params_derivs.m, ie use efficient sylvester equation method to compute analytical derivatives as in Ratto & Iskrev (2012)
options_mom_.initialize_estimated_parameters_with_the_prior_mode = 0; % needed by set_prior.m
options_mom_.figures = options_.figures; % needed by plot_priors.m
options_mom_.ramsey_policy = false; % needed by evaluate_steady_state
options_mom_.risky_steadystate = false; % needed by resol
options_mom_.jacobian_flag = true; % needed by dynare_solve
options_mom_.use_mh_covariance_matrix = false; % needed by posterior_sampler, get's overwritten by same option in options_mom_.posterior_sampler_options
% -------------------------------------------------------------------------
% CHECKS ON SETTINGS
% -------------------------------------------------------------------------
if strcmp(mom_method,'GMM') || strcmp(mom_method,'SMM')
if numel(options_mom_.nobs) > 1
error('method_of_moments: Recursive estimation is not supported. Please set an integer as ''nobs''!');
end
if numel(options_mom_.first_obs) > 1
error('method_of_moments: Recursive estimation is not supported. Please set an integer as ''first_obs''!');
end
end
if options_mom_.order < 1
error('method_of_moments: The order of the Taylor approximation cannot be 0!')
end
if options_mom_.order > 2
fprintf('Dynare will use ''k_order_solver'' as the order>2\n');
options_mom_.k_order_solver = true;
end
if strcmp(mom_method,'SMM')
if options_mom_.mom.simulation_multiple < 1
fprintf('The simulation horizon is shorter than the data. Dynare resets the simulation_multiple to 7.\n')
options_mom_.mom.simulation_multiple = 7;
end
end
if strcmp(mom_method,'GMM')
% require pruning with GMM at higher order
if options_mom_.order > 1 && ~options_mom_.pruning
fprintf('GMM at higher order only works with pruning, so we set pruning option to 1.\n');
options_mom_.pruning = true;
end
if options_mom_.order > 3
error('method_of_moments: Perturbation orders higher than 3 are not implemented for GMM estimation, try using SMM!');
end
end
if strcmp(mom_method,'IRF_MATCHING') && do_bayesian_estimation
if isfield(options_mom_,'mh_tune_jscale') && options_mom_.mh_tune_jscale.status && (options_mom_.mh_tune_jscale.maxiter<options_mom_.mh_tune_jscale.stepsize)
warning('method_of_moments: You specified mh_tune_jscale, but the maximum number of iterations is smaller than the step size. No update will take place.')
end
if options_mom_.load_results_after_load_mh
if ~exist([options_mom_.dirname filesep 'method_of_moments' filesep fname '_mom_results.mat'],'file')
fprintf('\nYou specified the ''load_results_after_load_mh'' option, but no ''%s_mom_results.mat'' file\n',fname);
fprintf('was found in the folder %s%smethod_of_moments.\n',options_mom_.dirname,filesep);
fprintf('Results will be recomputed and option ''load_results_after_load_mh'' is reset to false.\n');
options_mom_.load_results_after_load_mh = false;
end
end
if options_mom_.mh_replic>0 && options_mom_.mh_nblck<1
error('method_of_moments: Bayesian MCMC estimation cannot be conducted with ''mh_nblocks''=0!')
end
end
if options_mom_.mom.analytic_jacobian && ~strcmp(mom_method,'GMM')
options_mom_.mom.analytic_jacobian = false;
fprintf('\n''analytic_jacobian'' option will be dismissed as it only works with GMM.\n');
end
if strcmp(options_mom_.mom.mom_method,'IRF_MATCHING')
if any(cellfun(@(x) isnumeric(x) && any(x == 13), options_mom_.optimizer_vec))
error('method_of_moments: lsqnonlin (mode_compute=13) is not yet supported for IRF matching!');
end
end
\ No newline at end of file
matlab/+mom/display_comparison_moments.m → matlab/+mom/display_comparison_moments_irfs.m
View file @ 4a32acc2
function display_comparison_moments(M_, options_mom_, data_moments, model_moments)
% function display_comparison_moments(M_, options_mom_, data_moments, model_moments)
function display_comparison_moments_irfs(M_, options_mom_, data_moments, model_moments)
% display_comparison_moments_irfs(M_, options_mom_, data_moments, model_moments)
% -------------------------------------------------------------------------
% Displays and saves to disk the comparison of the data moments/IRFs and the model moments/IRFs
% -------------------------------------------------------------------------
% Displays and saves to disk the comparison of the data moments and the model moments
% =========================================================================
% INPUTS
% M_: [structure] model information
% options_mom_: [structure] method of moments options
......@@ -19,7 +19,8 @@ function display_comparison_moments(M_, options_mom_, data_moments, model_moment
% o dyn_latex_table
% o dyntable
% o cellofchararraymaxlength
% =========================================================================
% -------------------------------------------------------------------------
% Copyright © 2023 Dynare Team
%
% This file is part of Dynare.
......@@ -36,8 +37,28 @@ function display_comparison_moments(M_, options_mom_, data_moments, model_moment
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% =========================================================================
if strcmp(options_mom_.mom.mom_method,'IRF_MATCHING')
titl = upper('Comparison of matched data IRFs and model IRFs');
headers = {'IRF','Data','Model'};
idx = 1;
for jj = 1:size(M_.matched_irfs,1)
irf_varname = M_.matched_irfs{jj,1};
irf_shockname = M_.matched_irfs{jj,2};
% note that periods can span over multiple rows
IRF_PERIODS = [];
for kk = 1:size(M_.matched_irfs{jj,3},1)
irf_periods = M_.matched_irfs{jj,3}{kk,1};
IRF_PERIODS = [IRF_PERIODS; irf_periods(:)];
end
for hh = 1:length(IRF_PERIODS)
labels{idx,1} = sprintf('%s %s (%u)',irf_varname,irf_shockname,IRF_PERIODS(hh));
labels_TeX{idx,1} = sprintf('%s %s (%u)',M_.endo_names_tex{ismember(M_.endo_names,irf_varname)},M_.exo_names_tex{ismember(M_.exo_names,irf_shockname)},IRF_PERIODS(hh));
idx = idx+1;
end
end
else
titl = ['Comparison of matched data moments and model moments (',options_mom_.mom.mom_method,')'];
headers = {'Moment','Data','Model'};
for jm = 1:size(M_.matched_moments,1)
......@@ -67,6 +88,7 @@ for jm = 1:size(M_.matched_moments,1)
labels{jm,1} = lables_tmp;
labels_TeX{jm,1} = lables_tmp_tex;
end
end
data_mat = [data_moments model_moments];
dyntable(options_mom_, titl, headers, labels, data_mat, cellofchararraymaxlength(labels)+2, 10, 7);
if options_mom_.TeX
......
matlab/+mom/get_data_moments.m
View file @ 4a32acc2
function [dataMoments, m_data] = get_data_moments(data, oo_, matched_moments_, options_mom_)
% [dataMoments, m_data] = get_data_moments(data, oo_, matched_moments_, options_mom_)
% This function computes the user-selected empirical moments from data
% =========================================================================
function [data_moments, m_data] = get_data_moments(data, obs_var, inv_order_var, matched_moments_, options_mom_)
% [data_moments, m_data] = get_data_moments(data, obs_var, inv_order_var, matched_moments_, options_mom_)
% -------------------------------------------------------------------------
% Computes the user-selected empirical moments from data
% -------------------------------------------------------------------------
% INPUTS
% o data [T x varobs_nbr] data set
% o oo_: [structure] storage for results
% o obs_var: [integer] index of observables
% o inv_order_var: [integer] inverse decision rule order
% o matched_moments_: [structure] information about selected moments to match in estimation
% o options_mom_: [structure] information about all settings (specified by the user, preprocessor, and taken from global options_)
% -------------------------------------------------------------------------
% OUTPUTS
% o dataMoments [numMom x 1] mean of selected empirical moments
% o data_moments [numMom x 1] mean of selected empirical moments
% o m_data [T x numMom] selected empirical moments at each point in time
% -------------------------------------------------------------------------
% This function is called by
% o mom.run
% o mom.objective_function
% =========================================================================
% -------------------------------------------------------------------------
% Copyright © 2020-2023 Dynare Team
%
% This file is part of Dynare.
......@@ -32,24 +35,20 @@ function [dataMoments, m_data] = get_data_moments(data, oo_, matched_moments_, o
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------
% Author(s):
% o Willi Mutschler (willi@mutschler.eu)
% o Johannes Pfeifer (johannes.pfeifer@unibw.de)
% =========================================================================
% Initialization
T = size(data,1); % Number of observations (T)
dataMoments = NaN(options_mom_.mom.mom_nbr,1);
data_moments = NaN(options_mom_.mom.mom_nbr,1);
m_data = NaN(T,options_mom_.mom.mom_nbr);
% Product moment for each time period, i.e. each row t contains y_t1(l1)^p1*y_t2(l2)^p2*...
% note that here we already are able to treat leads and lags and any power product moments
for jm = 1:options_mom_.mom.mom_nbr
vars = oo_.dr.inv_order_var(matched_moments_{jm,1})';
vars = inv_order_var(matched_moments_{jm,1})';
leadlags = matched_moments_{jm,2}; % lags are negative numbers and leads are positive numbers
powers = matched_moments_{jm,3};
for jv = 1:length(vars)
jvar = (oo_.mom.obs_var == vars(jv));
jvar = (obs_var == vars(jv));
y = NaN(T,1); %Take care of T_eff instead of T for lags and NaN via mean with 'omitnan' option below
y( (1-min(leadlags(jv),0)) : (T-max(leadlags(jv),0)), 1) = data( (1+max(leadlags(jv),0)) : (T+min(leadlags(jv),0)), jvar).^powers(jv);
if jv==1
......@@ -59,11 +58,11 @@ for jm = 1:options_mom_.mom.mom_nbr
end
end
% We replace NaN (due to leads and lags and missing values) with the corresponding mean
if isoctave || matlab_ver_less_than('8.5')
dataMoments(jm,1) = nanmean(m_data_tmp);
if isoctave && octave_ver_less_than('8')
data_moments(jm,1) = nanmean(m_data_tmp);
else
dataMoments(jm,1) = mean(m_data_tmp,'omitnan');
data_moments(jm,1) = mean(m_data_tmp,'omitnan');
end
m_data_tmp(isnan(m_data_tmp)) = dataMoments(jm,1);
m_data_tmp(isnan(m_data_tmp)) = data_moments(jm,1);
m_data(:,jm) = m_data_tmp;
end
matlab/+mom/graph_comparison_irfs.m 0 → 100644
View file @ 4a32acc2
function graph_comparison_irfs(matched_irfs,irf_model_varobs,varobs_id,irf_horizon,relative_irf,endo_names,endo_names_tex,exo_names,exo_names_tex,dname,fname,graph_format,TeX,nodisplay,figures_textwidth)
% graph_comparison_irfs(matched_irfs,irf_model_varobs,varobs_id,irf_horizon,relative_irf,endo_names,endo_names_tex,exo_names,exo_names_tex,dname,fname,graph_format,TeX,nodisplay,figures_textwidth)
% -------------------------------------------------------------------------
% Plots and saves to disk the comparison of the selected data IRFs and corresponding model IRfs
% -------------------------------------------------------------------------
% INPUTS
% matched_irfs: [matrix] information on matched data IRFs
% irf_model_varobs: [matrix] model IRFs for observable variables
% varobs_id: [vector] index for observable variables in endo_names
% irf_horizon: [scalar] maximum horizon of IRFs
% relative_irf: [boolean] if true, plots normalized IRFs
% endo_names: [cell] names of endogenous variables
% endo_names_tex: [cell] names of endogenous variables in latex
% exo_names: [cell] names of exogenous variables
% exo_names_tex: [cell] names of exogenous variables in latex
% dname: [string] name of the directory where to save the graphs
% fname: [string] name of the mod file
% graph_format: [cell] format of the graphs
% TeX: [boolean] if true, uses latex for plots
% nodisplay: [boolean] if true, does not display the graphs
% figures_textwidth: [scalar] textwidth used in plots
% -------------------------------------------------------------------------
% OUTPUT
% No output, just displays and saves to disk the graphs
% -------------------------------------------------------------------------
% This function is called by
% o mom.run
% -------------------------------------------------------------------------
% This function calls
% o dyn_figure
% o dyn_saveas
% o remove_fractional_xticks
% o CheckPath
% o pltorg
% -------------------------------------------------------------------------
% Copyright © 2023 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
graph_directory_name = CheckPath('graphs',dname);
latex_directory_name = CheckPath('latex',dname);
if TeX && any(strcmp('eps',cellstr(graph_format)))
fid_TeX = fopen([latex_directory_name '/' fname '_irf_matching_plot.tex'],'w');
fprintf(fid_TeX,'%% TeX eps-loader file generated by mom.run.m (Dynare).\n');
fprintf(fid_TeX,['%% ' datestr(now,0) '\n']);
fprintf(fid_TeX,' \n');
end
unique_shock_entries = unique(matched_irfs(:, 2));
colDarkGrey = [0.3, 0.3, 0.3]; % dark grey
for jexo = unique_shock_entries' % loop over cell with shock names
unique_variables = unique(matched_irfs(ismember(matched_irfs(:, 2),jexo), 1));
[nbplt,nr,nc,lr,lc,nstar] = pltorg(length(unique_variables));
fig = 0;
for jvar = 1:length(unique_variables)
% get data points, note that periods and values can span over multiple rows
jj = ismember(matched_irfs(:,1), unique_variables(jvar)) & ismember(matched_irfs(:,2), jexo);
IRF_PERIODS = []; IRF_VALUES = [];
for kk = 1:size(matched_irfs{jj,3},1)
irf_periods = matched_irfs{jj,3}{kk,1};
irf_values = matched_irfs{jj,3}{kk,2};
if length(irf_values)==1
irf_values = repmat(irf_values,length(irf_periods),1);
end
IRF_PERIODS = [IRF_PERIODS; irf_periods(:)];
IRF_VALUES = [IRF_VALUES; irf_values(:)];
end
if jvar==1 || ~( (fig-1)*nstar<jvar && jvar<=fig*nstar )
fig = fig+1;
fig_irf = dyn_figure(nodisplay,'Name',['IRF matching shock to ' jexo{:} ' figure ' int2str(fig)]);
end
plt = jvar-(fig-1)*nstar;
if nbplt>1 && fig==nbplt
subplot(lr,lc,plt);
else
subplot(nr,nc,plt);
end
plt_data = plot(IRF_PERIODS,IRF_VALUES,'h', 'MarkerEdgeColor',colDarkGrey,'MarkerFaceColor',colDarkGrey,'MarkerSize',8);
hold on
plt_model = plot(1:irf_horizon, irf_model_varobs(:,varobs_id==find(ismember(endo_names,unique_variables(jvar))) , ismember(exo_names,jexo)),'-k','linewidth',2);
hold on
plot([1 irf_horizon],[0 0],'-r','linewidth',1);
hold off
xlim([1 irf_horizon]);
remove_fractional_xticks
if TeX
title(['$' endo_names_tex{ismember(endo_names,unique_variables(jvar))} '$'],'Interpreter','latex');
else
title(unique_variables{jvar},'Interpreter','none');
end
set(gca,'FontSize',12);
if (plt==nstar) || jvar==length(unique_variables)
% Adding a legend at the bottom
axes('Position',[0, 0, 1, 1],'Visible','off');
lgd = legend([plt_data,plt_model],{'Data', 'Model'}, 'Location', 'southeast','NumColumns',2,'FontSize',14);
if ~isoctave
lgd.Position = [0.37 0.01 lgd.Position(3) lgd.Position(4)];
end
dyn_saveas(fig_irf,[graph_directory_name filesep fname '_matched_irf_' jexo{:} int2str(fig)],nodisplay,graph_format);
if TeX && any(strcmp('eps',cellstr(graph_format)))
fprintf(fid_TeX,'\\begin{figure}[H]\n');
fprintf(fid_TeX,'\\centering \n');
fprintf(fid_TeX,'\\includegraphics[width=%2.2f\\textwidth]{%s_matched_irf_%s%s}\n',figures_textwidth*min(plt/nc,1),[graph_directory_name '/' fname],jexo{:},int2str(fig));
if relative_irf
fprintf(fid_TeX,'\\caption{Relative impulse response functions (orthogonalized shock to $%s$).}', jexo{:});
else
fprintf(fid_TeX,'\\caption{Impulse response functions (orthogonalized shock to $%s$).}', jexo{:});
end
fprintf(fid_TeX,'\\label{Fig:MatchedIRF:%s:%s}\n', jexo{:},int2str(fig));
fprintf(fid_TeX,'\\end{figure}\n');
fprintf(fid_TeX,' \n');
end
end
end
end
\ No newline at end of file
matlab/+mom/matched_irfs_blocks.m 0 → 100644
View file @ 4a32acc2
function [data_irfs, weight_mat, irf_index, max_irf_horizon] = matched_irfs_blocks(matched_irfs, matched_irfs_weight, varobs_id, obs_nbr, exo_nbr, endo_names, exo_names)
% [data_irfs, weight_mat, irf_index, max_irf_horizon] = matched_irfs_blocks(matched_irfs, matched_irfs_weight, varobs_id, obs_nbr, exo_nbr, endo_names, exo_names)
% -------------------------------------------------------------------------
% Checks and transforms matched_irfs and matched_irfs_weight blocks
% for further use in the estimation.
% -------------------------------------------------------------------------
% INPUTS
% matched_irfs: [cell array] original matched_irfs block
% matched_irfs_weight: [cell array] original matched_irfs_weight block
% varobs_id: [vector] index for observable variables in endo_names
% obs_nbr: [scalar] number of observable variables
% exo_nbr: [scalar] number of exogenous variables
% endo_names: [cell array] names of endogenous variables
% exo_names: [cell array] names of exogenous variables
% -------------------------------------------------------------------------
% OUTPUT
% data_irfs: [matrix] IRFs for VAROBS as declared in matched_irfs block
% weight_mat: [matrix] weighting matrix for IRFs as declared in matched_irfs_weight block
% irf_index: [vector] index for selecting specific IRFs from full IRF matrix of observables
% max_irf_horizon: [scalar] maximum IRF horizon as declared in matched_irfs block
% -------------------------------------------------------------------------
% This function is called by
% o mom.run
% -------------------------------------------------------------------------
% Copyright © 2023 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% note matched_irfs block:
% - each row in the cell contains a unique combination of var and varexo,
% however the third column in each row is a nested cell with information
% on periods, values and weights
% - periods, values and weights can span several rows with different lengths of entries
% - in some cases we need to duplicate values and/or weights
% - at the end we want to have everything vectorized and the same length
% get maximum IRF horizons
max_irf_horizon = [];
for jj = 1:size(matched_irfs,1)
max_irf_horizon = [max_irf_horizon; cell2mat(cellfun(@(c) c(:), matched_irfs{jj,3}(:,1), 'UniformOutput', false))];
end
max_irf_horizon = max(max_irf_horizon);
% create full matrix where 1st dimension are IRF periods, 2nd dimension are variables as declared in VAROBS, 3rd dimension are shocks
% idea: overwrite NaN values if they are declared in matched_irfs block; at the end the remaining NaN values will be removed
data_irfs = NaN(max_irf_horizon,obs_nbr,exo_nbr);
% create full empirical weighting matrix, identity matrix by default, i.e. all IRFs are equally important
% idea: first specify full matrix and then reduce it using only entries that are declared in matched_irfs block
weight_mat = speye(max_irf_horizon*obs_nbr*exo_nbr);
for jj = 1:size(matched_irfs,1)
id_var = find(ismember(endo_names,matched_irfs{jj,1}));
id_varobs = find(varobs_id==id_var,1);
id_shock = find(ismember(exo_names,matched_irfs{jj,2}));
if isempty(id_varobs)
skipline;
error('method_of_moments: You specified an IRF matching involving variable %s, but it is not declared as a varobs!',endo_names{id_var})
end
IRF_PERIODS = []; IRF_VALUES = []; IRF_WEIGHTS = [];
for kk = 1:size(matched_irfs{jj,3},1)
irf_periods = matched_irfs{jj,3}{kk,1};
if length(unique(irf_periods)) < length(irf_periods) % row-specific check for unique periods
error('method_of_moments: You specified an IRF matching involving variable %s and shock %s, but there were duplicate ''periods'' in the specification!',endo_names{id_var},exo_names{id_shock});
end
irf_values = matched_irfs{jj,3}{kk,2};
if length(irf_values)==1
irf_values = repmat(irf_values,length(irf_periods),1);
end
if length(irf_periods) ~= length(irf_values) % row-specific check for enough values
error('method_of_moments: You specified an IRF matching involving variable %s and shock %s, but the length of ''periods'' does not match the length of ''values''!',endo_names{id_var},exo_names{id_shock});
end
irf_weights = matched_irfs{jj,3}{kk,3};
if length(irf_weights)==1
irf_weights = repmat(irf_weights,length(irf_periods),1);
end
if length(irf_periods) ~= length(irf_weights) % row-specific check for enough weights
error('method_of_moments: You specified an IRF matching involving variable %s and shock %s, but the length of ''periods'' does not match the length of ''weights''!',endo_names{id_var},exo_names{id_shock});
end
IRF_PERIODS = [IRF_PERIODS; irf_periods(:)];
IRF_VALUES = [IRF_VALUES; irf_values(:)];
IRF_WEIGHTS = [IRF_WEIGHTS; irf_weights(:)];
end
if length(unique(irf_periods)) < length(irf_periods) % overall check for unique periods
error('method_of_moments: You specified an IRF matching involving variable %s and shock %s, but there were duplicate ''periods'' in the specification!',endo_names{id_var},exo_names{id_shock});
end
for hh = 1:length(IRF_PERIODS)
data_irfs(IRF_PERIODS(hh),id_varobs,id_shock) = IRF_VALUES(hh);
if IRF_WEIGHTS(hh) ~= 1
idweight_mat = sub2ind(size(data_irfs),IRF_PERIODS(hh),id_varobs,id_shock);
weight_mat(idweight_mat,idweight_mat) = IRF_WEIGHTS(hh);
end
end
end
% fine-tune weighting matrix using matched_irfs_weights
for jj = 1:size(matched_irfs_weight,1)
id_var1 = find(ismember(endo_names,matched_irfs_weight{jj,1}));
id_var2 = find(ismember(endo_names,matched_irfs_weight{jj,4}));
id_varobs1 = find(varobs_id==id_var1,1);
id_varobs2 = find(varobs_id==id_var2,1);
if isempty(id_varobs1)
skipline;
error('method_of_moments: You specified a weight for an IRF matching involving variable %s, but it is not a varobs!',endo_names{id_var1})
end
if isempty(id_varobs2)
skipline;
error('method_of_moments: You specified a weight for an IRF matching involving variable %s, but it is not a varobs!',endo_names{id_var2})
end
id_shock1 = find(ismember(exo_names,matched_irfs_weight{jj,3}));
id_shock2 = find(ismember(exo_names,matched_irfs_weight{jj,6}));
irf_periods1 = matched_irfs_weight{jj,2};
irf_periods2 = matched_irfs_weight{jj,5};
if length(irf_periods1) ~= length(irf_periods2)
error('method_of_moments: You specified a ''matched_irfs_weights'' entry for an IRF matching involving %s/%s and %s/%s,\n but the horizons do not have the same length!',endo_names{id_var1},exo_names{id_shock1},endo_names{id_var2},exo_names{id_shock2});
end
if max([irf_periods1(:);irf_periods2(:)]) > max_irf_horizon
error('method_of_moments: You specified a ''matched_irfs_weights'' entry for an IRF matching involving %s/%s and %s/%s,\n but the horizon is larger than the maximum one declared in the ''matched_irfs'' block!',endo_names{id_var1},exo_names{id_shock1},endo_names{id_var2},exo_names{id_shock2});
end
weight_mat_values = matched_irfs_weight{jj,7};
if length(weight_mat_values)==1 && length(irf_periods1)>1
weight_mat_values = repmat(weight_mat_values,length(irf_periods1),1);
end
if length(weight_mat_values) ~= length(irf_periods1)
error('method_of_moments: You specified a ''matched_irfs_weights'' entry for an IRF matching involving %s/%s and %s/%s,\n but the horizons do not match the length of ''weights''!',endo_names{id_var1},exo_names{id_shock1},endo_names{id_var2},exo_names{id_shock2});
end
for hh = 1:length(irf_periods1)
idweight_mat1 = sub2ind(size(data_irfs),irf_periods1(hh),id_varobs1,id_shock1);
idweight_mat2 = sub2ind(size(data_irfs),irf_periods2(hh),id_varobs2,id_shock2);
weight_mat(idweight_mat1,idweight_mat2) = weight_mat_values(hh);
weight_mat(idweight_mat2,idweight_mat1) = weight_mat_values(hh); % symmetry
end
end
% remove non-specified IRFs
irf_index = find(~isnan(data_irfs));
data_irfs = data_irfs(irf_index);
weight_mat = weight_mat(irf_index,irf_index);
\ No newline at end of file
matlab/+mom/matched_moments_block.m
View file @ 4a32acc2
function matched_moments = matched_moments_block(matched_moments, mom_method)
% function matched_moments = matched_moments_block(matched_moments, mom_method)
% matched_moments = matched_moments_block(matched_moments, mom_method)
% -------------------------------------------------------------------------
% Checks and transforms matched_moments block for further use in the estimation
% -------------------------------------------------------------------------
% Checks and transforms matched_moments bock for further use in the estimation
% =========================================================================
% INPUTS
% matched_moments: [cell array] original matched_moments block
% mom_method: [string] method of moments method (GMM or SMM)
......@@ -12,7 +12,8 @@ function matched_moments = matched_moments_block(matched_moments, mom_method)
% -------------------------------------------------------------------------
% This function is called by
% o mom.run
% =========================================================================
% -------------------------------------------------------------------------
% Copyright © 2023 Dynare Team
%
% This file is part of Dynare.
......@@ -29,7 +30,7 @@ function matched_moments = matched_moments_block(matched_moments, mom_method)
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% =========================================================================
matched_moments_orig = matched_moments;
% higher-order product moments not supported yet for GMM
......@@ -59,22 +60,22 @@ for jm = 1:size(matched_moments,1)
end
% find duplicate rows in cell array by making groups according to powers as we can then use cell2mat for the unique function
powers = cellfun(@sum,matched_moments(:,3))';
UniqueMomIdx = [];
unique_mom_idx = [];
for jpow = unique(powers)
idx1 = find(powers==jpow);
[~,idx2] = unique(cell2mat(matched_moments(idx1,:)),'rows');
UniqueMomIdx = [UniqueMomIdx idx1(idx2)];
unique_mom_idx = [unique_mom_idx idx1(idx2)];
end
% remove duplicate elements
DuplicateMoms = setdiff(1:size(matched_moments_orig,1),UniqueMomIdx);
if ~isempty(DuplicateMoms)
fprintf('Duplicate declared moments found and removed in ''matched_moments'' block in rows:\n %s.\n',num2str(DuplicateMoms))
duplicate_moms = setdiff(1:size(matched_moments_orig,1),unique_mom_idx);
if ~isempty(duplicate_moms)
fprintf('Duplicate declared moments found and removed in ''matched_moments'' block in rows:\n %s.\n',num2str(duplicate_moms));
fprintf('Dynare will continue with remaining moment conditions\n');
end
if strcmp(mom_method, 'SMM')
% for SMM: keep the original structure, but get rid of duplicate moments
matched_moments = matched_moments_orig(sort(UniqueMomIdx),:);
matched_moments = matched_moments_orig(sort(unique_mom_idx),:);
elseif strcmp(mom_method, 'GMM')
% for GMM we use the transformed matched_moments structure
matched_moments = matched_moments(sort(UniqueMomIdx),:);
matched_moments = matched_moments(sort(unique_mom_idx),:);
end
\ No newline at end of file
matlab/+mom/mode_compute_gmm_smm.m
View file @ 4a32acc2
function [xparam1, oo_, Woptflag] = mode_compute_gmm_smm(xparam0, objective_function, oo_, M_, options_mom_, estim_params_, bayestopt_, Bounds)
% function [xparam1, oo_, Woptflag] = mode_compute_gmm_smm(xparam0, objective_function, oo_, M_, options_mom_, estim_params_, bayestopt_, Bounds)
function [xparam1, weighting_info, mom_verbose] = mode_compute_gmm_smm(xparam0, objective_function, m_data, data_moments, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
% [xparam1, weighting_info, mom_verbose] = mode_compute_gmm_smm(xparam0, objective_function, m_data, data_moments, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
% -------------------------------------------------------------------------
% Iterated method of moments for GMM and SMM, computes the minimum of the
% objective function (distance between data moments and model moments)
% for a sequence of optimizers and GMM/SMM iterations with different
% weighting matrices.
% =========================================================================
% -------------------------------------------------------------------------
% INPUTS
% xparam0: [vector] vector of initialized parameters
% objective_function: [func handle] name of the objective function
% oo_: [structure] results
% M_: [structure] model information
% m_data: [matrix] selected data moments at each point in time
% data_moments: [vector] vector of data moments
% options_mom_: [structure] options
% M_: [structure] model information
% estim_params_: [structure] information on estimated parameters
% bayestopt_: [structure] information on priors
% Bounds: [structure] bounds for optimization
% BoundsInfo: [structure] bounds for optimization
% dr: [structure] reduced form model
% endo_steady_state: [vector] steady state for endogenous variables (initval)
% exo_steady_state: [vector] steady state for exogenous variables (initval)
% exo_det_steady_state: [vector] steady state for exogenous deterministic variables (initval)
% -------------------------------------------------------------------------
% OUTPUT
% xparam1: [vector] mode of objective function
% oo_: [structure] updated results
% Woptflag: [logical] true if optimal weighting matrix was computed
% weighting_info: [structure] information on weighting matrix
% mom_verbose: [structure] information on intermediate estimation results
% -------------------------------------------------------------------------
% This function is called by
% o mom.run
......@@ -29,7 +34,9 @@ function [xparam1, oo_, Woptflag] = mode_compute_gmm_smm(xparam0, objective_func
% o mom.display_estimation_results_table
% o dynare_minimize_objective
% o mom.objective_function
% =========================================================================
% o prior_dist_names
% -------------------------------------------------------------------------
% Copyright © 2023 Dynare Team
%
% This file is part of Dynare.
......@@ -46,15 +53,16 @@ function [xparam1, oo_, Woptflag] = mode_compute_gmm_smm(xparam0, objective_func
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% =========================================================================
mom_verbose = [];
if size(options_mom_.mom.weighting_matrix,1)>1 && ~(any(strcmpi('diagonal',options_mom_.mom.weighting_matrix)) || any(strcmpi('optimal',options_mom_.mom.weighting_matrix)))
fprintf('\nYou did not specify the use of an optimal or diagonal weighting matrix. There is no point in running an iterated method of moments.\n')
fprintf('\nYou did not specify the use of an optimal or diagonal weighting matrix. There is no point in running an iterated method of moments.\n');
end
for stage_iter = 1:size(options_mom_.mom.weighting_matrix,1)
fprintf('Estimation stage %u\n',stage_iter);
Woptflag = false;
weighting_info.Woptflag = false;
switch lower(options_mom_.mom.weighting_matrix{stage_iter})
case 'identity_matrix'
fprintf(' - identity weighting matrix\n');
......@@ -63,44 +71,44 @@ for stage_iter = 1:size(options_mom_.mom.weighting_matrix,1)
fprintf(' - diagonal of optimal weighting matrix (Bartlett kernel with %d lags)\n', options_mom_.mom.bartlett_kernel_lag);
if stage_iter == 1
fprintf(' and using data-moments as initial estimate of model-moments\n');
weighting_matrix = diag(diag( mom.optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.data_moments, options_mom_.mom.bartlett_kernel_lag) ));
weighting_matrix = diag(diag( mom.optimal_weighting_matrix(m_data, data_moments, options_mom_.mom.bartlett_kernel_lag) ));
else
fprintf(' and using previous stage estimate of model-moments\n');
weighting_matrix = diag(diag( mom.optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.model_moments, options_mom_.mom.bartlett_kernel_lag) ));
weighting_matrix = diag(diag( mom.optimal_weighting_matrix(m_data, model_moments, options_mom_.mom.bartlett_kernel_lag) ));
end
case 'optimal'
fprintf(' - optimal weighting matrix (Bartlett kernel with %d lags)\n', options_mom_.mom.bartlett_kernel_lag);
if stage_iter == 1
fprintf(' and using data-moments as initial estimate of model-moments\n');
weighting_matrix = mom.optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.data_moments, options_mom_.mom.bartlett_kernel_lag);
weighting_matrix = mom.optimal_weighting_matrix(m_data, data_moments, options_mom_.mom.bartlett_kernel_lag);
else
fprintf(' and using previous stage estimate of model-moments\n');
weighting_matrix = mom.optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.model_moments, options_mom_.mom.bartlett_kernel_lag);
Woptflag = true;
weighting_matrix = mom.optimal_weighting_matrix(m_data, model_moments, options_mom_.mom.bartlett_kernel_lag);
weighting_info.Woptflag = true;
end
otherwise % user specified matrix in file
fprintf(' - user-specified weighting matrix\n');
try
load(options_mom_.mom.weighting_matrix{stage_iter},'weighting_matrix')
catch
error(['method_of_moments: No matrix named ''weighting_matrix'' could be found in ',options_mom_.mom.weighting_matrix{stage_iter},'.mat !'])
error(['method_of_moments: No matrix named ''weighting_matrix'' could be found in ',options_mom_.mom.weighting_matrix{stage_iter},'.mat !']);
end
[nrow, ncol] = size(weighting_matrix);
if ~isequal(nrow,ncol) || ~isequal(nrow,length(oo_.mom.data_moments)) %check if square and right size
error(['method_of_moments: ''weighting_matrix'' must be square and have ',num2str(length(oo_.mom.data_moments)),' rows and columns!'])
if ~isequal(nrow,ncol) || ~isequal(nrow,length(data_moments)) %check if square and right size
error(['method_of_moments: ''weighting_matrix'' must be square and have ',num2str(length(data_moments)),' rows and columns!']);
end
end
try % check for positive definiteness of weighting_matrix
oo_.mom.Sw = chol(weighting_matrix);
weighting_info.Sw = chol(weighting_matrix);
catch
error('method_of_moments: Specified ''weighting_matrix'' is not positive definite. Check whether your model implies stochastic singularity!')
error('method_of_moments: Specified ''weighting_matrix'' is not positive definite. Check whether your model implies stochastic singularity!');
end
for optim_iter = 1:length(options_mom_.optimizer_vec)
options_mom_.current_optimizer = options_mom_.optimizer_vec{optim_iter};
if options_mom_.optimizer_vec{optim_iter} == 0
xparam1 = xparam0; % no minimization, evaluate objective at current values
fval = feval(objective_function, xparam1, Bounds, oo_, estim_params_, M_, options_mom_);
fval = feval(objective_function, xparam1, data_moments, weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
else
if options_mom_.optimizer_vec{optim_iter} == 13
options_mom_.mom.vector_output = true;
......@@ -112,18 +120,22 @@ for stage_iter = 1:size(options_mom_.mom.weighting_matrix,1)
else
options_mom_.mom.compute_derivs = false;
end
[xparam1, fval, exitflag] = dynare_minimize_objective(objective_function, xparam0, options_mom_.optimizer_vec{optim_iter}, options_mom_, [Bounds.lb Bounds.ub], bayestopt_.name, bayestopt_, [],...
Bounds, oo_, estim_params_, M_, options_mom_);
[xparam1, fval] = dynare_minimize_objective(objective_function, xparam0, options_mom_.optimizer_vec{optim_iter}, options_mom_, [BoundsInfo.lb BoundsInfo.ub], bayestopt_.name, bayestopt_, [],...
data_moments, weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
if options_mom_.mom.vector_output
fval = fval'*fval;
end
end
fprintf('\nStage %d Iteration %d: Value of minimized moment distance objective function: %12.10f.\n',stage_iter,optim_iter,fval)
fprintf('\nStage %d Iteration %d: Value of minimized moment distance objective function: %12.10f.\n',stage_iter,optim_iter,fval);
if options_mom_.mom.verbose
oo_.mom = display_estimation_results_table(xparam1,NaN(size(xparam1)),M_,options_mom_,estim_params_,bayestopt_,oo_.mom,prior_dist_names,sprintf('%s (STAGE %d ITERATION %d) VERBOSE',options_mom_.mom.mom_method,stage_iter,optim_iter),sprintf('verbose_%s_stage_%d_iter_%d',lower(options_mom_.mom.mom_method),stage_iter,optim_iter));
fprintf('\n''verbose'' option: ');
std_via_invhessian_xparam1_iter = NaN(size(xparam1));
tbl_title_iter = sprintf('FREQUENTIST %s (STAGE %d ITERATION %d) VERBOSE',options_mom_.mom.mom_method,stage_iter,optim_iter);
field_name_iter = sprintf('%s_stage_%d_iter_%d',lower(options_mom_.mom.mom_method),stage_iter,optim_iter);
mom_verbose.(field_name_iter) = display_estimation_results_table(xparam1,std_via_invhessian_xparam1_iter,M_,options_mom_,estim_params_,bayestopt_,[],prior_dist_names,tbl_title_iter,field_name_iter);
end
xparam0 = xparam1;
end
options_mom_.vector_output = false;
[~, ~, ~,~,~, oo_] = feval(objective_function, xparam1, Bounds, oo_, estim_params_, M_, options_mom_); % get oo_.mom.model_moments for iterated GMM/SMM to compute optimal weighting matrix
[~, ~, ~, ~, ~, ~, model_moments] = feval(objective_function, xparam1, data_moments, weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); % get model_moments for iterated GMM/SMM to compute optimal weighting matrix
end
matlab/+mom/mode_compute_irf_matching.m 0 → 100644
View file @ 4a32acc2
function [xparam1, hessian_xparam1, fval, mom_verbose] = mode_compute_irf_matching(xparam0, hessian_xparam0, objective_function, doBayesianEstimation, weighting_info, data_moments, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
% [xparam1, hessian_xparam1, fval, mom_verbose] = mode_compute_irf_matching(xparam0, hessian_xparam0, objective_function, doBayesianEstimation, weighting_info, data_moments, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
% -------------------------------------------------------------------------
% Computes the minimum of the objective function (distance between data IRFs
% and model IRFs) for a sequence of optimizers.
% Note that we call a "mode" the minimum of the objective function, i.e.
% the parameter vector that minimizes the distance between the IRFs
% computed from the model and the IRFs computed from the data.
% -------------------------------------------------------------------------
% INPUTS
% xparam0: [vector] initialized parameters
% hessian_xparam0: [matrix] initialized Hessian at xparam0
% objective_function: [func handle] name of the objective function
% doBayesianEstimation: [logical] true if Bayesian estimation
% weighting_info: [structure] information on weighting matrix
% data_moments: [vector] data moments
% options_mom_: [structure] options
% M_: [structure] model information
% estim_params_: [structure] information on estimated parameters
% bayestopt_: [structure] information on priors
% BoundsInfo: [structure] bounds for optimization
% dr: [structure] information reduced-form model
% endo_steady_state: [vector] steady state of endogenous variables (initval)
% exo_steady_state: [vector] steady state of exogenous variables (initval)
% exo_det_steady_state: [vector] steady state of deterministic exogenous variables (initval)
% -------------------------------------------------------------------------
% OUTPUT
% xparam1: [vector] mode of objective function
% hessian_xparam1: [matrix] Hessian at xparam1
% fval: [double] function value at mode
% mom_verbose: [structure] information on intermediate estimation results
% Also saves the computed mode and hessian to a file.
% -------------------------------------------------------------------------
% This function is called by
% o mom.run
% -------------------------------------------------------------------------
% This function calls
% o display_estimation_results_table
% o dynare_minimize_objective
% o hessian
% o mom.objective_function
% -------------------------------------------------------------------------
% Copyright © 2023 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
mom_verbose = [];
for optim_iter = 1:length(options_mom_.optimizer_vec)
options_mom_.current_optimizer = options_mom_.optimizer_vec{optim_iter};
if options_mom_.optimizer_vec{optim_iter}==0
% no minimization, evaluate objective at current values
xparam1 = xparam0;
hessian_xparam1 = hessian_xparam0;
fval = feval(objective_function, xparam1, data_moments, weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
else
[xparam1, fval, exitflag, hessian_xparam1, options_mom_, Scale, new_rat_hess_info] = dynare_minimize_objective(objective_function, xparam0, options_mom_.optimizer_vec{optim_iter}, options_mom_, [BoundsInfo.lb BoundsInfo.ub], bayestopt_.name, bayestopt_, hessian_xparam0,...
data_moments, weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
end
fprintf('\nMode Compute Iteration %d: Value of minimized moment distance objective function: %12.10f.\n',optim_iter,fval);
if options_mom_.mom.verbose
fprintf('\n''verbose'' option: ');
if options_mom_.cova_compute
if options_mom_.optimizer_vec{optim_iter}==0
hessian_xparam1_iter = hessian_xparam1;
else
fprintf('computing Hessian');
hessian_xparam1_iter = hessian(objective_function, xparam1, options_mom_.gstep,...
data_moments, weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
hessian_xparam1_iter = reshape(hessian_xparam1_iter, length(xparam1), length(xparam1));
end
hsd_iter = sqrt(diag(hessian_xparam1_iter));
invhessian_xparam1_iter = inv(hessian_xparam1_iter./(hsd_iter*hsd_iter'))./(hsd_iter*hsd_iter');
std_via_invhessian_xparam1_iter = sqrt(diag(invhessian_xparam1_iter));
else
std_via_invhessian_xparam1_iter = NaN(size(xparam1));
end
fprintf(' and displaying intermediate results.');
if doBayesianEstimation
tbl_title_iter = sprintf('BAYESIAN %s (OPTIM ITERATION %d) VERBOSE',strrep(options_mom_.mom.mom_method,'_',' '),optim_iter);
field_name_iter = sprintf('posterior_iter_%d',optim_iter);
else
tbl_title_iter = sprintf('FREQUENTIST %s (OPTIM ITERATION %d) VERBOSE',strrep(options_mom_.mom.mom_method,'_',' '),optim_iter);
field_name_iter = sprintf('iter_%d',optim_iter);
end
mom_verbose.(field_name_iter) = display_estimation_results_table(xparam1,std_via_invhessian_xparam1_iter,M_,options_mom_,estim_params_,bayestopt_,[],prior_dist_names,tbl_title_iter,field_name_iter);
end
xparam0 = xparam1;
hessian_xparam0 = hessian_xparam1;
end
if options_mom_.cova_compute
if options_mom_.mom.verbose
hessian_xparam1 = hessian_xparam1_iter;
else
fprintf('\nComputing Hessian at the mode.\n');
hessian_xparam1 = hessian(objective_function, xparam1, options_mom_.gstep,...
data_moments, weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
hessian_xparam1 = reshape(hessian_xparam1, length(xparam1), length(xparam1));
end
end
parameter_names = bayestopt_.name;
if options_mom_.cova_compute || options_mom_.mode_compute==5 || options_mom_.mode_compute==6
hh = hessian_xparam1;
save([M_.dname filesep 'method_of_moments' filesep M_.fname '_mode.mat'],'xparam1','hh','parameter_names','fval');
else
save([M_.dname filesep 'method_of_moments' filesep M_.fname '_mode.mat'],'xparam1','parameter_names','fval');
end
\ No newline at end of file
matlab/+mom/objective_function.m
View file @ 4a32acc2
function [fval, info, exit_flag, df, junkHessian, oo_, M_] = objective_function(xparam, Bounds, oo_, estim_params_, M_, options_mom_, bayestopt_)
% [fval, info, exit_flag, df, junk1, oo_, M_] = objective_function(xparam, Bounds, oo_, estim_params_, M_, options_mom_, bayestopt_)
function [fval, info, exit_flag, df, junk_hessian, Q, model_moments, model_moments_params_derivs, irf_model_varobs] = objective_function(xparam, data_moments, weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
% [fval, info, exit_flag, df, junk_hessian, Q, model_moments, model_moments_params_derivs, irf_model_varobs] = objective_function(xparam, data_moments, weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
% -------------------------------------------------------------------------
% This function evaluates the objective function for method of moments estimation
% =========================================================================
% INPUTS
% o xparam: [vector] current value of estimated parameters as returned by set_prior()
% o Bounds: [structure] containing parameter bounds
% o oo_: [structure] for results
% o estim_params_: [structure] describing the estimated_parameters
% o M_ [structure] describing the model
% o options_mom_: [structure] information about all settings (specified by the user, preprocessor, and taken from global options_)
% o bayestopt_: [structure] information about the prior
% This function evaluates the objective function for method of moments estimation,
% i.e. distance between model and data moments/IRFs, possibly augmented with priors.
% -------------------------------------------------------------------------
% INPUTS (same ones as in dsge_likelihood.m and dsge_var_likelihood.m)
% - xparam: [vector] current value of estimated parameters as returned by set_prior()
% - data_moments: [vector] data with moments/IRFs to match (corresponds to dataset_ in likelihood-based estimation)
% - weighting_info: [structure] storing information on weighting matrices (corresponds to dataset_info in likelihood-based estimation)
% - options_mom_: [structure] information about all settings (specified by the user, preprocessor, and taken from global options_)
% - M_ [structure] model information
% - estim_params_: [structure] information from estimated_params block
% - bayestopt_: [structure] information on the prior distributions
% - BoundsInfo: [structure] parameter bounds
% - dr: [structure] reduced form model
% - endo_steady_state: [vector] steady state value for endogenous variables (initval)
% - exo_steady_state: [vector] steady state value for exogenous variables (initval)
% - exo_det_steady_state: [vector] steady state value for exogenous deterministic variables (initval)
% -------------------------------------------------------------------------
% OUTPUTS
% o fval: [double] value of the quadratic form of the moment difference (except for lsqnonlin, where this is done implicitly)
% o info: [vector] information on error codes and penalties
% o exit_flag: [double] flag for exit status (0 if error, 1 if no error)
% o df: [matrix] analytical jacobian of the moment difference (wrt paramters), currently for GMM only
% o junkHessian: [matrix] empty matrix required for optimizer interface (Hessian would typically go here)
% o oo_: [structure] results with the following updated fields:
% - oo_.mom.model_moments: [vector] model moments
% - oo_.mom.Q: [double] value of the quadratic form of the moment difference
% - oo_.mom.model_moments_params_derivs: [matrix] analytical jacobian of the model moments wrt estimated parameters (currently for GMM only)
% o M_: [structure] updated model structure
% - fval: [double] value of the quadratic form of the moment difference (except for lsqnonlin, where this is done implicitly)
% - info: [vector] information on error codes and penalties
% - exit_flag: [double] flag for exit status (0 if error, 1 if no error)
% - df: [matrix] analytical Jacobian of the moment difference (wrt parameters), currently for GMM only
% - junk_hessian: [matrix] empty matrix required for optimizer interface (Hessian would typically go here)
% - Q: [double] value of the quadratic form of the moment difference
% - model_moments: [vector] model moments
% - model_moments_params_derivs: [matrix] analytical Jacobian of the model moments wrt estimated parameters (currently for GMM only)
% - irf_model_varobs: [matrix] model IRFs for observable variables (used for plotting matched IRfs in mom.run)
% -------------------------------------------------------------------------
% This function is called by
% o mom.run
......@@ -30,13 +35,18 @@ function [fval, info, exit_flag, df, junkHessian, oo_, M_] = objective_function(
% -------------------------------------------------------------------------
% This function calls
% o check_bounds_and_definiteness_estimation
% o get_lower_cholesky_covariance
% o get_perturbation_params_derivs
% o getIrfShocksIndx
% o irf
% o mom.get_data_moments
% o priordens
% o pruned_state_space_system
% o resol
% o set_all_parameters
% o simult_
% =========================================================================
% -------------------------------------------------------------------------
% Copyright © 2020-2023 Dynare Team
%
% This file is part of Dynare.
......@@ -53,28 +63,27 @@ function [fval, info, exit_flag, df, junkHessian, oo_, M_] = objective_function(
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% =========================================================================
%% TO DO
% check the info values and make use of meaningful penalties
% how do we do the penalty for the prior??
%------------------------------------------------------------------------------
% Initialization of the returned variables and others...
%------------------------------------------------------------------------------
junkHessian = [];
irf_model_varobs = [];
model_moments_params_derivs = [];
model_moments = [];
Q = [];
junk_hessian = [];
df = []; % required to be empty by e.g. newrat
if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM')
if options_mom_.mom.compute_derivs && options_mom_.mom.analytic_jacobian
if options_mom_.mom.vector_output == 1
if options_mom_.mom.penalized_estimator
df = nan(size(oo_.mom.data_moments,1)+length(xparam),length(xparam));
df = NaN(options_mom_.mom.mom_nbr+length(xparam),length(xparam));
else
df = nan(size(oo_.mom.data_moments,1),length(xparam));
df = NaN(options_mom_.mom.mom_nbr,length(xparam));
end
else
df = nan(length(xparam),1);
df = NaN(length(xparam),1);
end
end
end
......@@ -83,25 +92,23 @@ end
%--------------------------------------------------------------------------
% Get the structural parameters and define penalties
%--------------------------------------------------------------------------
% Ensure that xparam1 is a column vector; particleswarm.m requires this.
xparam = xparam(:);
M_ = set_all_parameters(xparam, estim_params_, M_);
[fval,info,exit_flag] = check_bounds_and_definiteness_estimation(xparam, M_, estim_params_, Bounds);
[fval,info,exit_flag] = check_bounds_and_definiteness_estimation(xparam, M_, estim_params_, BoundsInfo);
if info(1)
if options_mom_.mom.vector_output == 1 % lsqnonlin requires vector output
fval = ones(size(oo_.mom.data_moments,1),1)*options_mom_.huge_number;
fval = ones(options_mom_.mom.mom_nbr,1)*options_mom_.huge_number;
end
return
end
%--------------------------------------------------------------------------
% Call resol to compute steady state and model solution
% Compute steady state and model solution
%--------------------------------------------------------------------------
% Compute linear approximation around the deterministic steady state
[oo_.dr, info, M_.params] = resol(0, M_, options_mom_, oo_.dr ,oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
[dr, info, M_.params] = compute_decision_rules(M_,options_mom_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
% Return, with endogenous penalty when possible, if resol issues an error code
if info(1)
if info(1) == 3 || info(1) == 4 || info(1) == 5 || info(1)==6 ||info(1) == 19 ||...
......@@ -109,10 +116,14 @@ if info(1)
info(1) == 81 || info(1) == 84 || info(1) == 85 || info(1) == 86
% meaningful second entry of output that can be used
fval = Inf;
if ~isfinite(info(2))
info(4) = 0.1;
else
info(4) = info(2);
end
exit_flag = 0;
if options_mom_.mom.vector_output == 1 % lsqnonlin requires vector output
fval = ones(size(oo_.mom.data_moments,1),1)*options_mom_.huge_number;
fval = ones(options_mom_.mom.mom_nbr,1)*options_mom_.huge_number;
end
return
else
......@@ -120,7 +131,7 @@ if info(1)
info(4) = 0.1;
exit_flag = 0;
if options_mom_.mom.vector_output == 1 % lsqnonlin requires vector output
fval = ones(size(oo_.mom.data_moments,1),1)*options_mom_.huge_number;
fval = ones(options_mom_.mom.mom_nbr,1)*options_mom_.huge_number;
end
return
end
......@@ -140,49 +151,49 @@ if strcmp(options_mom_.mom.mom_method,'GMM')
if ~isempty(estim_params_.var_exo)
indpstderr = estim_params_.var_exo(:,1); % values correspond to varexo declaration order, row number corresponds to order in estimated_params
end
indpcorr=[]; % initialize matrix for corr paramters
indpcorr=[]; % initialize matrix for corr parameters
if ~isempty(estim_params_.corrx)
indpcorr = estim_params_.corrx(:,1:2); % values correspond to varexo declaration order, row number corresponds to order in estimated_params
end
if estim_params_.nvn || estim_params_.ncn % nvn is number of stderr parameters and ncn is number of corr parameters of measurement innovations as declared in estimated_params
error('Analytic computation of standard errrors does not (yet) support measurement errors.\nInstead, define them explicitly as varexo and provide measurement equations in the model definition.\nAlternatively, use numerical standard errors.')
error('Analytic computation of standard errrors does not (yet) support measurement errors.\nInstead, define them explicitly as varexo and provide measurement equations in the model definition.\nAlternatively, use numerical standard errors.');
end
modparam_nbr = estim_params_.np; % number of model parameters as declared in estimated_params
stderrparam_nbr = estim_params_.nvx; % number of stderr parameters
corrparam_nbr = estim_params_.ncx; % number of corr parameters
totparam_nbr = stderrparam_nbr+corrparam_nbr+modparam_nbr;
oo_.dr.derivs = get_perturbation_params_derivs(M_, options_mom_, estim_params_, oo_, indpmodel, indpstderr, indpcorr, 0); %analytic derivatives of perturbation matrices
oo_.mom.model_moments_params_derivs = NaN(options_mom_.mom.mom_nbr,totparam_nbr);
pruned_state_space = pruned_state_space_system(M_, options_mom_, oo_.dr, oo_.mom.obs_var, options_mom_.ar, 0, 1);
dr.derivs = identification.get_perturbation_params_derivs(M_, options_mom_, estim_params_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, indpmodel, indpstderr, indpcorr, 0); %analytic derivatives of perturbation matrices
model_moments_params_derivs = NaN(options_mom_.mom.mom_nbr,totparam_nbr);
pruned_state_space = pruned_SS.pruned_state_space_system(M_, options_mom_, dr, options_mom_.mom.obs_var, options_mom_.ar, 0, 1);
else
pruned_state_space = pruned_state_space_system(M_, options_mom_, oo_.dr, oo_.mom.obs_var, options_mom_.ar, 0, 0);
pruned_state_space = pruned_SS.pruned_state_space_system(M_, options_mom_, dr, options_mom_.mom.obs_var, options_mom_.ar, 0, 0);
end
oo_.mom.model_moments = NaN(options_mom_.mom.mom_nbr,1);
model_moments = NaN(options_mom_.mom.mom_nbr,1);
for jm = 1:size(M_.matched_moments,1)
% First moments
if ~options_mom_.prefilter && (sum(M_.matched_moments{jm,3}) == 1)
idx1 = (oo_.mom.obs_var == find(oo_.dr.order_var==M_.matched_moments{jm,1}) );
oo_.mom.model_moments(jm,1) = pruned_state_space.E_y(idx1);
idx1 = (options_mom_.mom.obs_var == find(dr.order_var==M_.matched_moments{jm,1}) );
model_moments(jm,1) = pruned_state_space.E_y(idx1);
if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
oo_.mom.model_moments_params_derivs(jm,:) = pruned_state_space.dE_y(idx1,:);
model_moments_params_derivs(jm,:) = pruned_state_space.dE_y(idx1,:);
end
end
% second moments
if (sum(M_.matched_moments{jm,3}) == 2)
idx1 = (oo_.mom.obs_var == find(oo_.dr.order_var==M_.matched_moments{jm,1}(1)) );
idx2 = (oo_.mom.obs_var == find(oo_.dr.order_var==M_.matched_moments{jm,1}(2)) );
idx1 = (options_mom_.mom.obs_var == find(dr.order_var==M_.matched_moments{jm,1}(1)) );
idx2 = (options_mom_.mom.obs_var == find(dr.order_var==M_.matched_moments{jm,1}(2)) );
if nnz(M_.matched_moments{jm,2}) == 0
% covariance
if options_mom_.prefilter
oo_.mom.model_moments(jm,1) = pruned_state_space.Var_y(idx1,idx2);
model_moments(jm,1) = pruned_state_space.Var_y(idx1,idx2);
if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
oo_.mom.model_moments_params_derivs(jm,:) = pruned_state_space.dVar_y(idx1,idx2,:);
model_moments_params_derivs(jm,:) = pruned_state_space.dVar_y(idx1,idx2,:);
end
else
oo_.mom.model_moments(jm,1) = pruned_state_space.Var_y(idx1,idx2) + pruned_state_space.E_y(idx1)*pruned_state_space.E_y(idx2)';
model_moments(jm,1) = pruned_state_space.Var_y(idx1,idx2) + pruned_state_space.E_y(idx1)*pruned_state_space.E_y(idx2)';
if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
for jp=1:totparam_nbr
oo_.mom.model_moments_params_derivs(jm,jp) = pruned_state_space.dVar_y(idx1,idx2,jp) + pruned_state_space.dE_y(idx1,jp)*pruned_state_space.E_y(idx2)' + pruned_state_space.E_y(idx1)*pruned_state_space.dE_y(idx2,jp)';
model_moments_params_derivs(jm,jp) = pruned_state_space.dVar_y(idx1,idx2,jp) + pruned_state_space.dE_y(idx1,jp)*pruned_state_space.E_y(idx2)' + pruned_state_space.E_y(idx1)*pruned_state_space.dE_y(idx2,jp)';
end
end
end
......@@ -190,15 +201,15 @@ if strcmp(options_mom_.mom.mom_method,'GMM')
% autocovariance
lag = -M_.matched_moments{jm,2}(2); %note that leads/lags in M_.matched_moments are transformed such that first entry is always 0 and the second is a lag
if options_mom_.prefilter
oo_.mom.model_moments(jm,1) = pruned_state_space.Var_yi(idx1,idx2,lag);
model_moments(jm,1) = pruned_state_space.Var_yi(idx1,idx2,lag);
if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
oo_.mom.model_moments_params_derivs(jm,:) = pruned_state_space.dVar_yi(idx1,idx2,lag,:);
model_moments_params_derivs(jm,:) = pruned_state_space.dVar_yi(idx1,idx2,lag,:);
end
else
oo_.mom.model_moments(jm,1) = pruned_state_space.Var_yi(idx1,idx2,lag) + pruned_state_space.E_y(idx1)*pruned_state_space.E_y(idx2)';
model_moments(jm,1) = pruned_state_space.Var_yi(idx1,idx2,lag) + pruned_state_space.E_y(idx1)*pruned_state_space.E_y(idx2)';
if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
for jp=1:totparam_nbr
oo_.mom.model_moments_params_derivs(jm,jp) = vec( pruned_state_space.dVar_yi(idx1,idx2,lag,jp) + pruned_state_space.dE_y(idx1,jp)*pruned_state_space.E_y(idx2)' + pruned_state_space.E_y(idx1)*pruned_state_space.dE_y(idx2,jp)');
model_moments_params_derivs(jm,jp) = vec( pruned_state_space.dVar_yi(idx1,idx2,lag,jp) + pruned_state_space.dE_y(idx1,jp)*pruned_state_space.E_y(idx2)' + pruned_state_space.E_y(idx1)*pruned_state_space.dE_y(idx2,jp)');
end
end
end
......@@ -218,26 +229,22 @@ if strcmp(options_mom_.mom.mom_method,'SMM')
scaled_shock_series = zeros(size(options_mom_.mom.shock_series)); % initialize
scaled_shock_series(:,i_exo_var) = options_mom_.mom.shock_series(:,i_exo_var)*chol_S; % set non-zero entries
% simulate series
y_sim = simult_(M_, options_mom_, oo_.dr.ys, oo_.dr, scaled_shock_series, options_mom_.order);
% provide meaningful penalty if data is nan or inf
y_sim = simult_(M_, options_mom_, dr.ys, dr, scaled_shock_series, options_mom_.order);
% provide meaningful penalty if data is NaN or Inf
if any(any(isnan(y_sim))) || any(any(isinf(y_sim)))
if options_mom_.mom.vector_output == 1 % lsqnonlin requires vector output
fval = Inf(size(oo_.mom.Sw,1),1);
else
fval = Inf;
end
info(1)=180;
info(4) = 0.1;
exit_flag = 0;
fval = Inf;
if options_mom_.mom.vector_output == 1 % lsqnonlin requires vector output
fval = ones(size(oo_.mom.data_moments,1),1)*options_mom_.huge_number;
fval = ones(options_mom_.mom.mom_nbr,1)*options_mom_.huge_number;
end
return
end
% remove burn-in and focus on observables (note that y_sim is in declaration order)
y_sim = y_sim(oo_.dr.order_var(oo_.mom.obs_var) , end-options_mom_.mom.long+1:end)';
y_sim = y_sim(dr.order_var(options_mom_.mom.obs_var) , end-options_mom_.mom.long+1:end)';
if ~all(diag(M_.H)==0)
i_ME = setdiff([1:size(M_.H,1)],find(diag(M_.H) == 0)); % find ME with 0 variance
i_ME = setdiff(1:size(M_.H,1),find(diag(M_.H) == 0)); % find ME with 0 variance
chol_S = chol(M_.H(i_ME,i_ME)); % decompose rest
shock_mat=zeros(size(options_mom_.mom.ME_shock_series)); % initialize
shock_mat(:,i_ME)=options_mom_.mom.ME_shock_series(:,i_ME)*chol_S;
......@@ -247,46 +254,144 @@ if strcmp(options_mom_.mom.mom_method,'SMM')
if options_mom_.prefilter
y_sim = bsxfun(@minus, y_sim, mean(y_sim,1));
end
oo_.mom.model_moments = mom.get_data_moments(y_sim, oo_, M_.matched_moments, options_mom_);
model_moments = mom.get_data_moments(y_sim, options_mom_.mom.obs_var, dr.inv_order_var, M_.matched_moments, options_mom_);
end
%------------------------------------------------------------------------------
% IRF_MATCHING using STOCH_SIMUL: Compute IRFs given model solution and Cholesky
% decomposition on shock covariance matrix; this resembles the core codes in
% stoch_simul
%------------------------------------------------------------------------------
if strcmp(options_mom_.mom.mom_method,'IRF_MATCHING') && strcmp(options_mom_.mom.simulation_method,'STOCH_SIMUL')
cs = get_lower_cholesky_covariance(M_.Sigma_e,options_mom_.add_tiny_number_to_cholesky);
irf_shocks_indx = getIrfShocksIndx(M_, options_mom_);
model_irf = NaN(options_mom_.irf,M_.endo_nbr,M_.exo_nbr);
for i = irf_shocks_indx
if options_mom_.order>1 && options_mom_.relative_irf % normalize shock to 0.01 before IRF generation for GIRFs; multiply with 100 later
y_irf = irf(M_, options_mom_, dr, cs(:,i)./cs(i,i)/100, options_mom_.irf, options_mom_.drop, options_mom_.replic, options_mom_.order);
else % for linear model, rescaling is done later
y_irf = irf(M_, options_mom_, dr, cs(:,i), options_mom_.irf, options_mom_.drop, options_mom_.replic, options_mom_.order);
end
if any(any(isnan(y_irf))) && ~options_mom_.pruning && ~(options_mom_.order==1)
info(1) = 181;
info(4) = 0.1;
fval = Inf;
exit_flag = 0;
if options_mom_.mom.vector_output == 1 % lsqnonlin requires vector output
fval = ones(options_mom_.mom.mom_nbr,1)*options_mom_.huge_number;
end
message = get_error_message(info,options_mom_);
fprintf('\n%s\n info = %d for shock %s.\n', message, info(1), M_.exo_names{i});
return
end
if options_mom_.relative_irf
if options_mom_.order==1 % multiply with 100 for backward compatibility
y_irf = 100*y_irf/cs(i,i);
end
end
model_irf(:,:,i) = transpose(y_irf);
end
% do transformations on model IRFs if irf_matching_file is provided
if ~isempty(options_mom_.mom.irf_matching_file.name)
[model_irf, check] = feval(str2func(options_mom_.mom.irf_matching_file.name), model_irf, M_, options_mom_, dr.ys);
if check
fval = Inf;
info(1) = 182;
info(4) = 0.1;
exit_flag = 0;
if options_mom_.mom.vector_output == 1 % lsqnonlin requires vector output
fval = ones(options_mom_.mom.mom_nbr,1)*options_mom_.huge_number;
end
return
end
end
irf_model_varobs = model_irf(:,options_mom_.varobs_id,:); % focus only on observables (this will be used later for plotting)
model_moments = irf_model_varobs(options_mom_.mom.irfIndex); % focus only on selected IRF periods
end
%--------------------------------------------------------------------------
% Compute quadratic target function
%--------------------------------------------------------------------------
moments_difference = oo_.mom.data_moments - oo_.mom.model_moments;
moments_difference = data_moments - model_moments;
if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM')
residuals = sqrt(options_mom_.mom.weighting_matrix_scaling_factor)*oo_.mom.Sw*moments_difference;
oo_.mom.Q = residuals'*residuals;
if strcmp(options_mom_.mom.mom_method,'IRF_MATCHING')
Q = transpose(moments_difference)*weighting_info.W*moments_difference;
% log-likelihood
lnlik = options_mom_.mom.mom_nbr/2*log(1/2/pi) - 1/2*weighting_info.Winv_logdet - 1/2*Q;
if isinf(lnlik)
fval = Inf; info(1) = 50; info(4) = 0.1; exit_flag = 0;
return
end
if isnan(lnlik)
fval = Inf; info(1) = 45; info(4) = 0.1; exit_flag = 0;
return
end
if imag(lnlik)~=0
fval = Inf; info(1) = 46; info(4) = 0.1; exit_flag = 0;
return
end
% add log prior if necessary
lnprior = priordens(xparam,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7,bayestopt_.p3,bayestopt_.p4);
if isinf(lnprior)
fval = Inf; info(1) = 40; info(4) = 0.1; exit_flag = 0;
return
end
if isnan(lnprior)
fval = Inf; info(1) = 47; info(4) = 0.1; exit_flag = 0;
return
end
if imag(lnprior)~=0
fval = Inf; info(1) = 48; info(4) = 0.1; exit_flag = 0;
return
end
% compute the posterior kernel
fval = - (lnlik + lnprior);
elseif strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM')
residuals = sqrt(options_mom_.mom.weighting_matrix_scaling_factor)*weighting_info.Sw*moments_difference;
Q = residuals'*residuals;
if isinf(Q)
fval = Inf; info(1) = 177; info(4) = 0.1; exit_flag = 0;
return
end
if isnan(Q)
fval = Inf; info(1) = 178; info(4) = 0.1; exit_flag = 0;
return
end
if imag(Q)~=0
fval = Inf; info(1) = 179; info(4) = 0.1; exit_flag = 0;
return
end
if options_mom_.mom.vector_output == 1 % lsqnonlin requires vector output
fval = residuals;
if options_mom_.mom.penalized_estimator
fval=[fval;(xparam-oo_.mom.prior.mean)./sqrt(diag(oo_.mom.prior.variance))];
fval=[fval;(xparam-bayestopt_.p1)./bayestopt_.p2];
end
else
fval = oo_.mom.Q;
fval = Q;
if options_mom_.mom.penalized_estimator
fval=fval+(xparam-oo_.mom.prior.mean)'/oo_.mom.prior.variance*(xparam-oo_.mom.prior.mean);
fval=fval+(xparam-bayestopt_.p1)'/(diag(bayestopt_.p2.^2))*(xparam-bayestopt_.p1);
end
end
if options_mom_.mom.compute_derivs && options_mom_.mom.analytic_jacobian
if options_mom_.mom.penalized_estimator
dxparam1 = eye(length(xparam));
dxparam = eye(length(xparam));
end
for jp=1:length(xparam)
dmoments_difference = - oo_.mom.model_moments_params_derivs(:,jp);
dresiduals = sqrt(options_mom_.mom.weighting_matrix_scaling_factor)*oo_.mom.Sw*dmoments_difference;
dmoments_difference = - model_moments_params_derivs(:,jp);
dresiduals = sqrt(options_mom_.mom.weighting_matrix_scaling_factor)*weighting_info.Sw*dmoments_difference;
if options_mom_.mom.vector_output == 1 % lsqnonlin requires vector output
if options_mom_.mom.penalized_estimator
df(:,jp)=[dresiduals;dxparam1(:,jp)./sqrt(diag(oo_.mom.prior.variance))];
df(:,jp)=[dresiduals;dxparam(:,jp)./bayestopt_.p2];
else
df(:,jp) = dresiduals;
end
else
df(jp,1) = dresiduals'*residuals + residuals'*dresiduals;
if options_mom_.mom.penalized_estimator
df(jp,1)=df(jp,1)+(dxparam1(:,jp))'/oo_.mom.prior.variance*(xparam-oo_.mom.prior.mean)+(xparam-oo_.mom.prior.mean)'/oo_.mom.prior.variance*(dxparam1(:,jp));
df(jp,1)=df(jp,1)+(dxparam(:,jp))'/(diag(bayestopt_.p2.^2))*(xparam-bayestopt_.p1)+(xparam-bayestopt_.p1)'/(diag(bayestopt_.p2.^2))*(dxparam(:,jp));
end
end
end
......@@ -295,4 +400,3 @@ end
end % main function end
\ No newline at end of file
matlab/+mom/optimal_weighting_matrix.m
View file @ 4a32acc2
......@@ -2,9 +2,10 @@ function W_opt = optimal_weighting_matrix(m_data, moments, q_lag)
% W_opt = optimal_weighting_matrix(m_data, moments, q_lag)
% -------------------------------------------------------------------------
% This function computes the optimal weigthing matrix by a Bartlett kernel with maximum lag q_lag
% Adapted from replication codes of
% o Andreasen, Fernández-Villaverde, Rubio-Ramírez (2018): "The Pruned State-Space System for Non-Linear DSGE Models: Theory and Empirical Applications", Review of Economic Studies, 85(1):1-49.
% =========================================================================
% Adapted from replication codes of Andreasen, Fernández-Villaverde, Rubio-Ramírez (2018):
% "The Pruned State-Space System for Non-Linear DSGE Models: Theory and Empirical Applications",
% Review of Economic Studies, 85(1):1-49.
% -------------------------------------------------------------------------
% INPUTS
% o m_data [T x numMom] selected data moments at each point in time
% o moments [numMom x 1] selected estimated moments (either data_moments or estimated model_moments)
......@@ -17,9 +18,10 @@ function W_opt = optimal_weighting_matrix(m_data, moments, q_lag)
% o mom.run.m
% -------------------------------------------------------------------------
% This function calls:
% o CorrMatrix (embedded)
% =========================================================================
% Copyright © 2020-2021 Dynare Team
% o corr_matrix (embedded)
% -------------------------------------------------------------------------
% Copyright © 2020-2023 Dynare Team
%
% This file is part of Dynare.
%
......@@ -35,46 +37,42 @@ function W_opt = optimal_weighting_matrix(m_data, moments, q_lag)
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------
% Author(s):
% o Willi Mutschler (willi@mutschler.eu)
% o Johannes Pfeifer (johannes.pfeifer@unibw.de)
% =========================================================================
% Initialize
% initialize
[T,num_Mom] = size(m_data); % note that in m_data NaN values (due to leads or lags in matched_moments and missing data) were replaced by the mean
% center around moments (could be either data_moments or model_moments)
h_Func = m_data - repmat(moments',T,1);
h_func = m_data - repmat(moments',T,1);
% The required correlation matrices
GAMA_array = zeros(num_Mom,num_Mom,q_lag);
GAMA0 = Corr_Matrix(h_Func,T,num_Mom,0);
% the required correlation matrices
gamma_array = zeros(num_Mom,num_Mom,q_lag);
gamma0 = corr_matrix(h_func,T,num_Mom,0);
if q_lag > 0
for ii=1:q_lag
GAMA_array(:,:,ii) = Corr_Matrix(h_Func,T,num_Mom,ii);
gamma_array(:,:,ii) = corr_matrix(h_func,T,num_Mom,ii);
end
end
% The estimate of S
S = GAMA0;
% the estimate of S
S = gamma0;
if q_lag > 0
for ii=1:q_lag
S = S + (1-ii/(q_lag+1))*(GAMA_array(:,:,ii) + GAMA_array(:,:,ii)');
S = S + (1-ii/(q_lag+1))*(gamma_array(:,:,ii) + gamma_array(:,:,ii)');
end
end
% The estimate of W
% the estimate of W
W_opt = S\eye(size(S,1));
W_opt=(W_opt+W_opt')/2; %assure symmetry
end
W_opt = (W_opt+W_opt')/2; % ensure symmetry
end % main function end
% The correlation matrix
function GAMA_corr = Corr_Matrix(h_Func,T,num_Mom,v)
GAMA_corr = zeros(num_Mom,num_Mom);
function gamma_corr = corr_matrix(h_func,T,num_Mom,v)
gamma_corr = zeros(num_Mom,num_Mom);
for t = 1+v:T
GAMA_corr = GAMA_corr + h_Func(t-v,:)'*h_Func(t,:);
end
GAMA_corr = GAMA_corr/T;
gamma_corr = gamma_corr + h_func(t-v,:)'*h_func(t,:);
end
gamma_corr = gamma_corr/T;
end % corr_matrix end
\ No newline at end of file