From 2e73856f5a4381c760ab1f3a8d9c423aa0702d46 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Thu, 7 Dec 2023 21:15:18 +0100 Subject: [PATCH 1/3] GSA and identification: move files to namespace --- license.txt | 44 +++++----- matlab/{gsa => +gsa}/Morris_Measure_Groups.m | 0 matlab/{gsa => +gsa}/Sampling_Function_2.m | 0 matlab/{gsa/myboxplot.m => +gsa/boxplot.m} | 4 +- matlab/{gsa => +gsa}/cumplot.m | 0 .../log_trans_.m => +gsa/log_transform.m} | 8 +- matlab/{gsa => +gsa}/map_calibration.m | 16 ++-- .../map_identification.m} | 44 +++++----- .../monte_carlo_filtering.m} | 50 +++++------ .../monte_carlo_filtering_analysis.m} | 18 ++-- .../monte_carlo_moments.m} | 10 +-- .../prior_draw_gsa.m => +gsa/prior_draw.m} | 3 +- matlab/{gsa => +gsa}/priorcdf.m | 0 .../reduced_form_mapping.m} | 40 ++++----- .../reduced_form_screening.m} | 16 ++-- matlab/{dynare_sensitivity.m => +gsa/run.m} | 18 ++-- matlab/{gsa => +gsa}/scatter_analysis.m | 4 +- matlab/{gsa => +gsa}/scatter_mcf.m | 4 +- matlab/{gsa => +gsa}/scatter_plots.m | 2 +- matlab/{gsa => +gsa}/set_shocks_param.m | 0 .../{gsa/gsa_skewness.m => +gsa/skewness.m} | 4 +- matlab/{gsa/smirnov.m => +gsa/smirnov_test.m} | 6 +- .../stab_map_.m => +gsa/stability_mapping.m} | 22 ++--- .../stability_mapping_bivariate.m} | 4 +- .../stability_mapping_univariate.m} | 16 ++-- .../stand_.m => +gsa/standardize_columns.m} | 4 +- matlab/{gsa => +gsa}/tcrit.m | 0 matlab/{gsa => +gsa}/teff.m | 0 matlab/{gsa => +gsa}/th_moments.m | 0 .../analysis.m} | 56 ++++++------- .../bruteforce.m} | 4 +- .../checks.m} | 12 +-- .../checks_via_subsets.m} | 8 +- matlab/{ => +identification}/cosn.m | 2 +- .../display.m} | 8 +- .../get_jacobians.m} | 12 +-- .../numerical_objective.m} | 4 +- .../plot.m} | 26 +++--- .../run.m} | 82 +++++++++---------- .../simulated_moment_uncertainty.m | 0 matlab/commutation.m | 2 +- matlab/duplication.m | 2 +- matlab/fjaco.m | 6 +- matlab/get_minimal_state_representation.m | 2 +- matlab/get_perturbation_params_derivs.m | 2 +- matlab/list_of_functions_to_be_cleared.m | 2 +- matlab/pruned_state_space_system.m | 4 +- matlab/set_all_parameters.m | 2 +- preprocessor | 2 +- 49 files changed, 288 insertions(+), 287 deletions(-) rename matlab/{gsa => +gsa}/Morris_Measure_Groups.m (100%) rename matlab/{gsa => +gsa}/Sampling_Function_2.m (100%) rename matlab/{gsa/myboxplot.m => +gsa/boxplot.m} (97%) rename matlab/{gsa => +gsa}/cumplot.m (100%) rename matlab/{gsa/log_trans_.m => +gsa/log_transform.m} (94%) rename matlab/{gsa => +gsa}/map_calibration.m (97%) rename matlab/{gsa/map_ident_.m => +gsa/map_identification.m} (90%) rename matlab/{gsa/filt_mc_.m => +gsa/monte_carlo_filtering.m} (95%) rename matlab/{gsa/mcf_analysis.m => +gsa/monte_carlo_filtering_analysis.m} (75%) rename matlab/{gsa/mc_moments.m => +gsa/monte_carlo_moments.m} (84%) rename matlab/{gsa/prior_draw_gsa.m => +gsa/prior_draw.m} (96%) rename matlab/{gsa => +gsa}/priorcdf.m (100%) rename matlab/{gsa/redform_map.m => +gsa/reduced_form_mapping.m} (95%) rename matlab/{gsa/redform_screen.m => +gsa/reduced_form_screening.m} (92%) rename matlab/{dynare_sensitivity.m => +gsa/run.m} (95%) rename matlab/{gsa => +gsa}/scatter_analysis.m (89%) rename matlab/{gsa => +gsa}/scatter_mcf.m (98%) rename matlab/{gsa => +gsa}/scatter_plots.m (99%) rename matlab/{gsa => +gsa}/set_shocks_param.m (100%) rename matlab/{gsa/gsa_skewness.m => +gsa/skewness.m} (94%) rename matlab/{gsa/smirnov.m => +gsa/smirnov_test.m} (94%) rename matlab/{gsa/stab_map_.m => +gsa/stability_mapping.m} (95%) rename matlab/{gsa/stab_map_2.m => +gsa/stability_mapping_bivariate.m} (96%) rename matlab/{gsa/stab_map_1.m => +gsa/stability_mapping_univariate.m} (88%) rename matlab/{gsa/stand_.m => +gsa/standardize_columns.m} (93%) rename matlab/{gsa => +gsa}/tcrit.m (100%) rename matlab/{gsa => +gsa}/teff.m (100%) rename matlab/{gsa => +gsa}/th_moments.m (100%) rename matlab/{identification_analysis.m => +identification/analysis.m} (94%) rename matlab/{ident_bruteforce.m => +identification/bruteforce.m} (98%) rename matlab/{identification_checks.m => +identification/checks.m} (95%) rename matlab/{identification_checks_via_subsets.m => +identification/checks_via_subsets.m} (98%) rename matlab/{ => +identification}/cosn.m (98%) rename matlab/{disp_identification.m => +identification/display.m} (98%) rename matlab/{get_identification_jacobians.m => +identification/get_jacobians.m} (97%) rename matlab/{identification_numerical_objective.m => +identification/numerical_objective.m} (97%) rename matlab/{plot_identification.m => +identification/plot.m} (95%) rename matlab/{dynare_identification.m => +identification/run.m} (94%) rename matlab/{ => +identification}/simulated_moment_uncertainty.m (100%) diff --git a/license.txt b/license.txt index d6f1336af..c1bee92b0 100644 --- a/license.txt +++ b/license.txt @@ -161,33 +161,33 @@ Comment: The author gave authorization to redistribute Journal of Multivariate Analysis, 2008, vol. 99, issue 3, pages 542-554. -Files: matlab/gsa/Morris_Measure_Groups.m - matlab/gsa/Sampling_Function_2.m +Files: matlab/+gsa/Morris_Measure_Groups.m + matlab/+gsa/Sampling_Function_2.m Copyright: 2005 European Commission - 2012-2017 Dynare Team + 2012-2013 Dynare Team License: GPL-3+ Comment: Written by Jessica Cariboni and Francesca Campolongo Joint Research Centre, The European Commission, -Files: matlab/gsa/cumplot.m - matlab/gsa/filt_mc_.m - matlab/gsa/gsa_skewness.m - matlab/gsa/log_trans_.m - matlab/gsa/map_calibration.m - matlab/gsa/map_ident_.m - matlab/gsa/mcf_analysis.m - matlab/gsa/myboxplot.m - matlab/gsa/prior_draw_gsa.m - matlab/gsa/redform_map.m - matlab/gsa/redform_screen.m - matlab/gsa/scatter_mcf.m - matlab/gsa/smirnov.m - matlab/gsa/stab_map_.m - matlab/gsa/stab_map_1.m - matlab/gsa/stab_map_2.m - matlab/gsa/stand_.m - matlab/gsa/tcrit.m - matlab/gsa/teff.m +Files: matlab/+gsa/cumplot.m + matlab/+gsa/monte_carlo_filtering.m + matlab/+gsa/skewness.m + matlab/+gsa/log_trans_.m + matlab/+gsa/map_calibration.m + matlab/+gsa/map_identification.m + matlab/+gsa/monte_carlo_filtering_analysis.m + matlab/+gsa/boxplot.m + matlab/+gsa/prior_draw.m + matlab/+gsa/reduced_form_mapping.m + matlab/+gsa/reduced_form_screening.m + matlab/+gsa/scatter_mcf.m + matlab/+gsa/smirnov_test.m + matlab/+gsa/stability_mapping.m + matlab/+gsa/stability_mapping_univariate.m + matlab/+gsa/stability_mapping_bivariate.m + matlab/+gsa/standardize_columns.m + matlab/+gsa/tcrit.m + matlab/+gsa/teff.m Copyright: 2011-2018 European Commission 2011-2023 Dynare Team License: GPL-3+ diff --git a/matlab/gsa/Morris_Measure_Groups.m b/matlab/+gsa/Morris_Measure_Groups.m similarity index 100% rename from matlab/gsa/Morris_Measure_Groups.m rename to matlab/+gsa/Morris_Measure_Groups.m diff --git a/matlab/gsa/Sampling_Function_2.m b/matlab/+gsa/Sampling_Function_2.m similarity index 100% rename from matlab/gsa/Sampling_Function_2.m rename to matlab/+gsa/Sampling_Function_2.m diff --git a/matlab/gsa/myboxplot.m b/matlab/+gsa/boxplot.m similarity index 97% rename from matlab/gsa/myboxplot.m rename to matlab/+gsa/boxplot.m index 4d6cf60d1..f893b7a81 100644 --- a/matlab/gsa/myboxplot.m +++ b/matlab/+gsa/boxplot.m @@ -1,5 +1,5 @@ -function sout = myboxplot (data,notched,symbol,vertical,maxwhisker) -% sout = myboxplot (data,notched,symbol,vertical,maxwhisker) +function sout = boxplot (data,notched,symbol,vertical,maxwhisker) +% sout = boxplot (data,notched,symbol,vertical,maxwhisker) % Creates a box plot % Copyright © 2010-2023 Dynare Team diff --git a/matlab/gsa/cumplot.m b/matlab/+gsa/cumplot.m similarity index 100% rename from matlab/gsa/cumplot.m rename to matlab/+gsa/cumplot.m diff --git a/matlab/gsa/log_trans_.m b/matlab/+gsa/log_transform.m similarity index 94% rename from matlab/gsa/log_trans_.m rename to matlab/+gsa/log_transform.m index 3dedb694e..852ddb187 100644 --- a/matlab/gsa/log_trans_.m +++ b/matlab/+gsa/log_transform.m @@ -1,5 +1,5 @@ -function [yy, xdir, isig, lam]=log_trans_(y0,xdir0,isig,lam) -% [yy, xdir, isig, lam]=log_trans_(y0,xdir0,isig,lam) +function [yy, xdir, isig, lam]=log_transform(y0,xdir0,isig,lam) +% [yy, xdir, isig, lam]=log_transform(y0,xdir0,isig,lam) % Conduct automatic log transformation lam(yy/isig+lam) % Inputs: % - y0 [double] series to transform @@ -56,10 +56,10 @@ end if nargin==1 xdir0=''; end -f=@(lam,y)gsa_skewness(log(y+lam)); +f=@(lam,y)gsa.skewness(log(y+lam)); isig=1; if ~(max(y0)<0 || min(y0)>0) - if gsa_skewness(y0)<0 + if gsa.skewness(y0)<0 isig=-1; y0=-y0; end diff --git a/matlab/gsa/map_calibration.m b/matlab/+gsa/map_calibration.m similarity index 97% rename from matlab/gsa/map_calibration.m rename to matlab/+gsa/map_calibration.m index 44703dc05..aa68e0900 100644 --- a/matlab/gsa/map_calibration.m +++ b/matlab/+gsa/map_calibration.m @@ -229,7 +229,7 @@ if ~isempty(indx_irf) if ~options_.nograph && length(time_matrix{plot_indx(ij)})==1 set(0,'currentfigure',h1), subplot(nrow,ncol, plot_indx(ij)), - hc = cumplot(mat_irf{ij}(:,ik)); + hc = gsa.cumplot(mat_irf{ij}(:,ik)); a=axis; delete(hc); x1val=max(endo_prior_restrictions.irf{ij,4}(1),a(1)); @@ -237,7 +237,7 @@ if ~isempty(indx_irf) hp = patch([x1val x2val x2val x1val],a([3 3 4 4]),'b'); hold all, set(hp,'FaceColor', [0.7 0.8 1]) - hc = cumplot(mat_irf{ij}(:,ik)); + hc = gsa.cumplot(mat_irf{ij}(:,ik)); set(hc,'color','k','linewidth',2) hold off, % hold off, @@ -259,7 +259,7 @@ if ~isempty(indx_irf) end options_mcf.title = atitle0; if ~isempty(indx1) && ~isempty(indx2) - mcf_analysis(xmat(:,nshock+1:end), indx1, indx2, options_mcf, M_, options_, bayestopt_, estim_params_); + gsa.monte_carlo_filtering_analysis(xmat(:,nshock+1:end), indx1, indx2, options_mcf, M_, options_, bayestopt_, estim_params_); end end for ij=1:nbr_irf_couples @@ -316,7 +316,7 @@ if ~isempty(indx_irf) options_mcf.title = atitle0; if ~isempty(indx1) && ~isempty(indx2) - mcf_analysis(xmat(:,nshock+1:end), indx1, indx2, options_mcf, M_, options_, bayestopt_, estim_params_); + gsa.monte_carlo_filtering_analysis(xmat(:,nshock+1:end), indx1, indx2, options_mcf, M_, options_, bayestopt_, estim_params_); end end end @@ -434,7 +434,7 @@ if ~isempty(indx_moment) if ~options_.nograph && length(time_matrix{plot_indx(ij)})==1 set(0,'currentfigure',h2); subplot(nrow,ncol,plot_indx(ij)), - hc = cumplot(mat_moment{ij}(:,ik)); + hc = gsa.cumplot(mat_moment{ij}(:,ik)); a=axis; delete(hc), % hist(mat_moment{ij}), x1val=max(endo_prior_restrictions.moment{ij,4}(1),a(1)); @@ -442,7 +442,7 @@ if ~isempty(indx_moment) hp = patch([x1val x2val x2val x1val],a([3 3 4 4]),'b'); set(hp,'FaceColor', [0.7 0.8 1]) hold all - hc = cumplot(mat_moment{ij}(:,ik)); + hc = gsa.cumplot(mat_moment{ij}(:,ik)); set(hc,'color','k','linewidth',2) hold off title([endo_prior_restrictions.moment{ij,1},' vs ',endo_prior_restrictions.moment{ij,2},'(',leg,')'],'interpreter','none'), @@ -463,7 +463,7 @@ if ~isempty(indx_moment) end options_mcf.title = atitle0; if ~isempty(indx1) && ~isempty(indx2) - mcf_analysis(xmat, indx1, indx2, options_mcf, M_, options_, bayestopt_, estim_params_); + gsa.monte_carlo_filtering_analysis(xmat, indx1, indx2, options_mcf, M_, options_, bayestopt_, estim_params_); end end for ij=1:nbr_moment_couples @@ -520,7 +520,7 @@ if ~isempty(indx_moment) end options_mcf.title = atitle0; if ~isempty(indx1) && ~isempty(indx2) - mcf_analysis(xmat, indx1, indx2, options_mcf, M_, options_, bayestopt_, estim_params_); + gsa.monte_carlo_filtering_analysis(xmat, indx1, indx2, options_mcf, M_, options_, bayestopt_, estim_params_); end end end diff --git a/matlab/gsa/map_ident_.m b/matlab/+gsa/map_identification.m similarity index 90% rename from matlab/gsa/map_ident_.m rename to matlab/+gsa/map_identification.m index 2b7194fa2..3f9f86b9b 100644 --- a/matlab/gsa/map_ident_.m +++ b/matlab/+gsa/map_identification.m @@ -1,5 +1,5 @@ -function map_ident_(OutputDirectoryName,opt_gsa,M_,oo_,options_,estim_params_,bayestopt_) -% map_ident_(OutputDirectoryName,opt_gsa,M_,oo_,options_,estim_params_,bayestopt_) +function map_identification(OutputDirectoryName,opt_gsa,M_,oo_,options_,estim_params_,bayestopt_) +% map_identification(OutputDirectoryName,opt_gsa,M_,oo_,options_,estim_params_,bayestopt_) % Inputs % - OutputDirectoryName [string] name of the output directory % - opt_gsa [structure] GSA options structure @@ -58,16 +58,16 @@ fname_ = M_.fname; if opt_gsa.load_ident_files==0 mss = yys(bayestopt_.mfys,:); - mss = teff(mss(:,istable),Nsam,istable); - yys = teff(yys(dr.order_var,istable),Nsam,istable); + mss = gsa.teff(mss(:,istable),Nsam,istable); + yys = gsa.teff(yys(dr.order_var,istable),Nsam,istable); if exist('T','var') - [vdec, cc, ac] = mc_moments(T, lpmatx, dr, M_, options_, estim_params_); + [vdec, cc, ac] = gsa.monte_carlo_moments(T, lpmatx, dr, M_, options_, estim_params_); else return end if opt_gsa.morris==2 - pdraws = dynare_identification(M_,oo_,options_,bayestopt_,estim_params_,options_.options_ident,[lpmatx lpmat(istable,:)]); + pdraws = identification.run(M_,oo_,options_,bayestopt_,estim_params_,options_.options_ident,[lpmatx lpmat(istable,:)]); if ~isempty(pdraws) && max(max(abs(pdraws-[lpmatx lpmat(istable,:)])))==0 disp(['Sample check OK. Largest difference: ', num2str(max(max(abs(pdraws-[lpmatx lpmat(istable,:)]))))]), clear pdraws; @@ -84,7 +84,7 @@ if opt_gsa.load_ident_files==0 end iplo=iplo+1; subplot(2,3,iplo) - myboxplot(squeeze(vdec(:,j,:))',[],'.',[],10) + gsa.boxplot(squeeze(vdec(:,j,:))',[],'.',[],10) set(gca,'xticklabel',' ','fontsize',10,'xtick',1:size(options_.varobs,1)) set(gca,'xlim',[0.5 size(options_.varobs,1)+0.5]) set(gca,'ylim',[-2 102]) @@ -105,11 +105,11 @@ if opt_gsa.load_ident_files==0 end end for j=1:size(cc,1) - cc(j,j,:)=stand_(squeeze(log(cc(j,j,:))))./2; + cc(j,j,:)=gsa.standardize_columns(squeeze(log(cc(j,j,:))))./2; end - [vdec, ~, ir_vdec, ic_vdec] = teff(vdec,Nsam,istable); - [cc, ~, ir_cc, ic_cc] = teff(cc,Nsam,istable); - [ac, ~, ir_ac, ic_ac] = teff(ac,Nsam,istable); + [vdec, ~, ir_vdec, ic_vdec] = gsa.teff(vdec,Nsam,istable); + [cc, ~, ir_cc, ic_cc] = gsa.teff(cc,Nsam,istable); + [ac, ~, ir_ac, ic_ac] = gsa.teff(ac,Nsam,istable); nc1= size(T,2); endo_nbr = M_.endo_nbr; @@ -123,7 +123,7 @@ if opt_gsa.load_ident_files==0 [Aa,Bb] = kalman_transition_matrix(dr,iv,ic); A = zeros(size(Aa,1),size(Aa,2)+size(Aa,1),length(istable)); if ~isempty(lpmatx) - M_=set_shocks_param(M_,estim_params_,lpmatx(1,:)); + M_=gsa.set_shocks_param(M_,estim_params_,lpmatx(1,:)); end A(:,:,1)=[Aa, triu(Bb*M_.Sigma_e*Bb')]; for j=2:length(istable) @@ -131,14 +131,14 @@ if opt_gsa.load_ident_files==0 dr.ghu = T(:, (nc1-M_.exo_nbr+1):end, j); [Aa,Bb] = kalman_transition_matrix(dr, iv, ic); if ~isempty(lpmatx) - M_=set_shocks_param(M_,estim_params_,lpmatx(j,:)); + M_=gsa.set_shocks_param(M_,estim_params_,lpmatx(j,:)); end A(:,:,j)=[Aa, triu(Bb*M_.Sigma_e*Bb')]; end clear T clear lpmatx - [yt, j0]=teff(A,Nsam,istable); + [yt, j0]=gsa.teff(A,Nsam,istable); yt = [yys yt]; if opt_gsa.morris==2 clear TAU A @@ -155,7 +155,7 @@ if opt_gsa.morris==1 if opt_gsa.load_ident_files==0 SAMorris=NaN(npT,3,size(vdec,2)); for i=1:size(vdec,2) - [~, SAMorris(:,:,i)] = Morris_Measure_Groups(npT, [lpmat0 lpmat], vdec(:,i),nliv); + [~, SAMorris(:,:,i)] = gsa.Morris_Measure_Groups(npT, [lpmat0 lpmat], vdec(:,i),nliv); end SAvdec = squeeze(SAMorris(:,1,:))'; save([OutputDirectoryName,'/',fname_,'_morris_IDE.mat'],'SAvdec','vdec','ir_vdec','ic_vdec') @@ -164,7 +164,7 @@ if opt_gsa.morris==1 end hh_fig = dyn_figure(options_.nodisplay,'name','Screening identification: variance decomposition'); - myboxplot(SAvdec,[],'.',[],10) + gsa.boxplot(SAvdec,[],'.',[],10) set(gca,'xticklabel',' ','fontsize',10,'xtick',1:npT) set(gca,'xlim',[0.5 npT+0.5]) ydum = get(gca,'ylim'); @@ -190,7 +190,7 @@ if opt_gsa.morris==1 ccac = [mss cc ac]; SAMorris=NaN(npT,3,size(ccac,2)); for i=1:size(ccac,2) - [~, SAMorris(:,:,i)] = Morris_Measure_Groups(npT, [lpmat0 lpmat], [ccac(:,i)],nliv); + [~, SAMorris(:,:,i)] = gsa.Morris_Measure_Groups(npT, [lpmat0 lpmat], [ccac(:,i)],nliv); end SAcc = squeeze(SAMorris(:,1,:))'; SAcc = SAcc./(max(SAcc,[],2)*ones(1,npT)); @@ -202,7 +202,7 @@ if opt_gsa.morris==1 end hh_fig=dyn_figure(options_.nodisplay,'name','Screening identification: theoretical moments'); - myboxplot(SAcc,[],'.',[],10) + gsa.boxplot(SAcc,[],'.',[],10) set(gca,'xticklabel',' ','fontsize',10,'xtick',1:npT) set(gca,'xlim',[0.5 npT+0.5]) set(gca,'ylim',[0 1]) @@ -223,7 +223,7 @@ if opt_gsa.morris==1 if opt_gsa.load_ident_files==0 SAMorris=NaN(npT,3,j0); for j=1:j0 - [~, SAMorris(:,:,j)] = Morris_Measure_Groups(npT, [lpmat0 lpmat], yt(:,j),nliv); + [~, SAMorris(:,:,j)] = gsa.Morris_Measure_Groups(npT, [lpmat0 lpmat], yt(:,j),nliv); end SAM = squeeze(SAMorris(1:end,1,:)); @@ -249,7 +249,7 @@ if opt_gsa.morris==1 load([OutputDirectoryName,'/',fname_,'_morris_IDE'],'SAnorm') end hh_fig=dyn_figure(options_.nodisplay,'name','Screening identification: model'); - myboxplot(SAnorm',[],'.',[],10) + gsa.boxplot(SAnorm',[],'.',[],10) set(gca,'xticklabel',' ','fontsize',10,'xtick',1:npT) set(gca,'xlim',[0.5 npT+0.5]) set(gca,'ylim',[0 1]) @@ -297,7 +297,7 @@ else % main effects analysis catch EET=[]; end - ccac = stand_([mss cc ac]); + ccac = gsa.standardize_columns([mss cc ac]); [pcc, dd] = eig(cov(ccac(istable,:))); [latent, isort] = sort(-diag(dd)); latent = -latent; @@ -314,7 +314,7 @@ else % main effects analysis if itrans==0 y0 = ccac(istable,j); elseif itrans==1 - y0 = log_trans_(ccac(istable,j)); + y0 = gsa.log_transform(ccac(istable,j)); else y0 = trank(ccac(istable,j)); end diff --git a/matlab/gsa/filt_mc_.m b/matlab/+gsa/monte_carlo_filtering.m similarity index 95% rename from matlab/gsa/filt_mc_.m rename to matlab/+gsa/monte_carlo_filtering.m index b59f6026c..69b100b4d 100644 --- a/matlab/gsa/filt_mc_.m +++ b/matlab/+gsa/monte_carlo_filtering.m @@ -1,5 +1,5 @@ -function [rmse_MC, ixx] = filt_mc_(OutDir,options_gsa_,dataset_,dataset_info,M_,oo_,options_,bayestopt_,estim_params_) -% [rmse_MC, ixx] = filt_mc_(OutDir,options_gsa_,dataset_,dataset_info,M_,oo_,options_,bayestopt_,estim_params_ +function [rmse_MC, ixx] = monte_carlo_filtering(OutDir,options_gsa_,dataset_,dataset_info,M_,oo_,options_,bayestopt_,estim_params_) +% [rmse_MC, ixx] = monte_carlo_filtering(OutDir,options_gsa_,dataset_,dataset_info,M_,oo_,options_,bayestopt_,estim_params_ % Inputs: % - OutputDirectoryName [string] name of the output directory % - options_gsa_ [structure] GSA options @@ -288,7 +288,7 @@ options_scatter.OutputDirectoryName = OutDir; options_scatter.amcf_name = asname; options_scatter.amcf_title = atitle; options_scatter.title = tmp_title; -scatter_analysis(r2_MC, x,options_scatter, options_); +gsa.scatter_analysis(r2_MC, x,options_scatter, options_); % end of visual scatter analysis if ~options_.opt_gsa.ppost && options_.opt_gsa.lik_only @@ -320,7 +320,7 @@ if ~options_.opt_gsa.ppost && options_.opt_gsa.lik_only options_mcf.nobeha_title_latex = 'worse posterior kernel'; end - mcf_analysis(x, ipost(1:nfilt), ipost(nfilt+1:end), options_mcf, M_, options_, bayestopt_, estim_params_); + gsa.monte_carlo_filtering_analysis(x, ipost(1:nfilt), ipost(nfilt+1:end), options_mcf, M_, options_, bayestopt_, estim_params_); if options_.opt_gsa.pprior anam = 'rmse_prior_lik'; atitle = 'RMSE prior: Log Likelihood Kernel'; @@ -338,7 +338,7 @@ if ~options_.opt_gsa.ppost && options_.opt_gsa.lik_only options_mcf.nobeha_title_latex = 'worse likelihood'; end - mcf_analysis(x, ilik(1:nfilt), ilik(nfilt+1:end), options_mcf, M_, options_, bayestopt_, estim_params_); + gsa.monte_carlo_filtering_analysis(x, ilik(1:nfilt), ilik(nfilt+1:end), options_mcf, M_, options_, bayestopt_, estim_params_); else if options_.opt_gsa.ppost @@ -367,9 +367,9 @@ else SS = zeros(npar+nshock, length(vvarvecm)); for j = 1:npar+nshock for i = 1:length(vvarvecm) - [~, P] = smirnov(x(ixx(nfilt0(i)+1:end,i),j),x(ixx(1:nfilt0(i),i),j), alpha); - [H1] = smirnov(x(ixx(nfilt0(i)+1:end,i),j),x(ixx(1:nfilt0(i),i),j),alpha,1); - [H2] = smirnov(x(ixx(nfilt0(i)+1:end,i),j),x(ixx(1:nfilt0(i),i),j),alpha,-1); + [~, P] = gsa.smirnov_test(x(ixx(nfilt0(i)+1:end,i),j),x(ixx(1:nfilt0(i),i),j), alpha); + [H1] = gsa.smirnov_test(x(ixx(nfilt0(i)+1:end,i),j),x(ixx(1:nfilt0(i),i),j),alpha,1); + [H2] = gsa.smirnov_test(x(ixx(nfilt0(i)+1:end,i),j),x(ixx(1:nfilt0(i),i),j),alpha,-1); if H1==0 && H2==0 SS(j,i)=1; elseif H1==0 @@ -382,7 +382,7 @@ else for i = 1:length(vvarvecm) for l = 1:length(vvarvecm) if l~=i && PP(j,i)10 && length(inobeha)>10 if options_.TeX - indcorr1 = stab_map_2(lpmat(ibeha,:),alpha2, pvalue_corr, M_, options_, bayestopt_, estim_params_, beha_title, beha_title_latex); - indcorr2 = stab_map_2(lpmat(inobeha,:),alpha2, pvalue_corr, M_, options_, bayestopt_, estim_params_, nobeha_title, nobeha_title_latex); + indcorr1 = gsa.stability_mapping_bivariate(lpmat(ibeha,:),alpha2, pvalue_corr, M_, options_, bayestopt_, estim_params_, beha_title, beha_title_latex); + indcorr2 = gsa.stability_mapping_bivariate(lpmat(inobeha,:),alpha2, pvalue_corr, M_, options_, bayestopt_, estim_params_, nobeha_title, nobeha_title_latex); else - indcorr1 = stab_map_2(lpmat(ibeha,:),alpha2, pvalue_corr, M_, options_, bayestopt_, estim_params_, beha_title); - indcorr2 = stab_map_2(lpmat(inobeha,:),alpha2, pvalue_corr, M_, options_, bayestopt_, estim_params_, nobeha_title); + indcorr1 = gsa.stability_mapping_bivariate(lpmat(ibeha,:),alpha2, pvalue_corr, M_, options_, bayestopt_, estim_params_, beha_title); + indcorr2 = gsa.stability_mapping_bivariate(lpmat(inobeha,:),alpha2, pvalue_corr, M_, options_, bayestopt_, estim_params_, nobeha_title); end indcorr = union(indcorr1(:), indcorr2(:)); indcorr = indcorr(~ismember(indcorr(:),indmcf)); @@ -104,11 +104,11 @@ if ~isempty(indmcf) && ~options_.nograph xx=xparam1(indmcf); end if options_.TeX - scatter_mcf(lpmat(ibeha,indmcf),lpmat(inobeha,indmcf), param_names_tex(indmcf), ... + gsa.scatter_mcf(lpmat(ibeha,indmcf),lpmat(inobeha,indmcf), param_names_tex(indmcf), ... '.', [fname_,'_',amcf_name], OutputDirectoryName, amcf_title,xx, options_, ... beha_title, nobeha_title, beha_title_latex, nobeha_title_latex) else - scatter_mcf(lpmat(ibeha,indmcf),lpmat(inobeha,indmcf), param_names_tex(indmcf), ... + gsa.scatter_mcf(lpmat(ibeha,indmcf),lpmat(inobeha,indmcf), param_names_tex(indmcf), ... '.', [fname_,'_',amcf_name], OutputDirectoryName, amcf_title,xx, options_, ... beha_title, nobeha_title) end diff --git a/matlab/gsa/mc_moments.m b/matlab/+gsa/monte_carlo_moments.m similarity index 84% rename from matlab/gsa/mc_moments.m rename to matlab/+gsa/monte_carlo_moments.m index 31554da12..510e28e94 100644 --- a/matlab/gsa/mc_moments.m +++ b/matlab/+gsa/monte_carlo_moments.m @@ -1,5 +1,5 @@ -function [vdec, cc, ac] = mc_moments(mm, ss, dr, M_, options_, estim_params_) -% [vdec, cc, ac] = mc_moments(mm, ss, dr, M_, options_,estim_params_) +function [vdec, cc, ac] = monte_carlo_moments(mm, ss, dr, M_, options_, estim_params_) +% [vdec, cc, ac] = monte_carlo_moments(mm, ss, dr, M_, options_,estim_params_) % Conduct Monte Carlo simulation of second moments for GSA % Inputs: % - dr [structure] decision rules @@ -32,7 +32,7 @@ function [vdec, cc, ac] = mc_moments(mm, ss, dr, M_, options_, estim_params_) [~, nc1, nsam] = size(mm); nobs=length(options_.varobs); -disp('mc_moments: Computing theoretical moments ...') +disp('monte_carlo_moments: Computing theoretical moments ...') h = dyn_waitbar(0,'Theoretical moments ...'); vdec = zeros(nobs,M_.exo_nbr,nsam); cc = zeros(nobs,nobs,nsam); @@ -42,9 +42,9 @@ for j=1:nsam dr.ghx = mm(:, 1:(nc1-M_.exo_nbr),j); dr.ghu = mm(:, (nc1-M_.exo_nbr+1):end, j); if ~isempty(ss) - M_=set_shocks_param(M_,estim_params_,ss(j,:)); + M_=gsa.set_shocks_param(M_,estim_params_,ss(j,:)); end - [vdec(:,:,j), corr, autocorr] = th_moments(dr,options_,M_); + [vdec(:,:,j), corr, autocorr] = gsa.th_moments(dr,options_,M_); cc(:,:,j)=triu(corr); dum=NaN(nobs,nobs*options_.ar); for i=1:options_.ar diff --git a/matlab/gsa/prior_draw_gsa.m b/matlab/+gsa/prior_draw.m similarity index 96% rename from matlab/gsa/prior_draw_gsa.m rename to matlab/+gsa/prior_draw.m index 58731ec0a..c3b8f8d9d 100644 --- a/matlab/gsa/prior_draw_gsa.m +++ b/matlab/+gsa/prior_draw.m @@ -1,4 +1,5 @@ -function pdraw = prior_draw_gsa(M_,bayestopt_,options_,estim_params_,init,rdraw) +function pdraw = prior_draw(M_,bayestopt_,options_,estim_params_,init,rdraw) +% pdraw = prior_draw(M_,bayestopt_,options_,estim_params_,init,rdraw) % Draws from the prior distributions for use with Sensitivity Toolbox for DYNARE % % INPUTS diff --git a/matlab/gsa/priorcdf.m b/matlab/+gsa/priorcdf.m similarity index 100% rename from matlab/gsa/priorcdf.m rename to matlab/+gsa/priorcdf.m diff --git a/matlab/gsa/redform_map.m b/matlab/+gsa/reduced_form_mapping.m similarity index 95% rename from matlab/gsa/redform_map.m rename to matlab/+gsa/reduced_form_mapping.m index 8475fd7e1..de575f4cb 100644 --- a/matlab/gsa/redform_map.m +++ b/matlab/+gsa/reduced_form_mapping.m @@ -1,5 +1,5 @@ -function redform_map(dirname,options_gsa_,M_,estim_params_,options_,bayestopt_,oo_) -% redform_map(dirname,options_gsa_,M_,estim_params_,options_,bayestopt_,oo_) +function reduced_form_mapping(dirname,options_gsa_,M_,estim_params_,options_,bayestopt_,oo_) +% reduced_form_mapping(dirname,options_gsa_,M_,estim_params_,options_,bayestopt_,oo_) % Inputs: % - dirname [string] name of the output directory % - options_gsa_ [structure] GSA options_ @@ -85,7 +85,7 @@ options_mcf.fname_ = M_.fname; options_mcf.OutputDirectoryName = adir; if ~exist('T','var') - stab_map_(dirname,options_gsa_,M_,oo_,options_,bayestopt_,estim_params_); + gsa.stability_mapping(dirname,options_gsa_,M_,oo_,options_,bayestopt_,estim_params_); if pprior load([dirname,filesep,M_.fname,'_prior'],'T'); else @@ -182,14 +182,14 @@ for j = 1:length(anamendo) end if ~options_.nograph hf=dyn_figure(options_.nodisplay,'name',['Reduced Form Mapping (Monte Carlo Filtering): ',namendo,' vs ', namexo]); - hc = cumplot(y0); + hc = gsa.cumplot(y0); a=axis; delete(hc); x1val=max(threshold(1),a(1)); x2val=min(threshold(2),a(2)); hp = patch([x1val x2val x2val x1val],a([3 3 4 4]),'b'); set(hp,'FaceColor', [0.7 0.8 1]) hold all, - hc = cumplot(y0); + hc = gsa.cumplot(y0); set(hc,'color','k','linewidth',2) hold off, if options_.TeX @@ -218,7 +218,7 @@ for j = 1:length(anamendo) options_mcf.OutputDirectoryName = xdir; if ~isempty(iy) && ~isempty(iyc) fprintf(['%4.1f%% of the ',type,' support matches ',atitle0,'\n'],length(iy)/length(y0)*100) - icheck = mcf_analysis(x0, iy, iyc, options_mcf, M_, options_, bayestopt_, estim_params_); + icheck = gsa.monte_carlo_filtering_analysis(x0, iy, iyc, options_mcf, M_, options_, bayestopt_, estim_params_); lpmat=x0(iy,:); if nshocks @@ -349,14 +349,14 @@ for j = 1:length(anamendo) end if ~options_.nograph hf=dyn_figure(options_.nodisplay,'name',['Reduced Form Mapping (Monte Carlo Filtering): ',namendo,' vs lagged ', namlagendo]); - hc = cumplot(y0); + hc = gsa.cumplot(y0); a=axis; delete(hc); x1val=max(threshold(1),a(1)); x2val=min(threshold(2),a(2)); hp = patch([x1val x2val x2val x1val],a([3 3 4 4]),'b'); set(hp,'FaceColor', [0.7 0.8 1]) hold all, - hc = cumplot(y0); + hc = gsa.cumplot(y0); set(hc,'color','k','linewidth',2) hold off if options_.TeX @@ -387,7 +387,7 @@ for j = 1:length(anamendo) if ~isempty(iy) && ~isempty(iyc) fprintf(['%4.1f%% of the ',type,' support matches ',atitle0,'\n'],length(iy)/length(y0)*100) - icheck = mcf_analysis(x0, iy, iyc, options_mcf, M_, options_, bayestopt_, estim_params_); + icheck = gsa.monte_carlo_filtering_analysis(x0, iy, iyc, options_mcf, M_, options_, bayestopt_, estim_params_); lpmat=x0(iy,:); if nshocks @@ -476,9 +476,9 @@ end if isempty(threshold) && ~options_.nograph hh_fig=dyn_figure(options_.nodisplay,'name','Reduced Form GSA'); if ilog==0 - myboxplot(si',[],'.',[],10) + gsa.boxplot(si',[],'.',[],10) else - myboxplot(silog',[],'.',[],10) + gsa.boxplot(silog',[],'.',[],10) end xlabel(' ') set(gca,'xticklabel',' ','fontsize',10,'xtick',1:np) @@ -513,7 +513,7 @@ if options_map.prior_range x0(:,j)=(x0(:,j)-pd(j,3))./(pd(j,4)-pd(j,3)); end else - x0=priorcdf(x0,pshape, pd(:,1), pd(:,2), pd(:,3), pd(:,4)); + x0=gsa.priorcdf(x0,pshape, pd(:,1), pd(:,2), pd(:,3), pd(:,4)); end if ilog @@ -549,7 +549,7 @@ if iload==0 ipred = setdiff(1:nrun,ifit); if ilog - [~, ~, isig, lam] = log_trans_(y0(iest)); + [~, ~, isig, lam] = gsa.log_transform(y0(iest)); y1 = log(y0*isig+lam); end if ~options_.nograph @@ -571,9 +571,9 @@ if iload==0 title(options_map.title,'interpreter','none') subplot(222) if ilog - hc = cumplot(y1); + hc = gsa.cumplot(y1); else - hc = cumplot(y0); + hc = gsa.cumplot(y0); end set(hc,'color','k','linewidth',2) title([options_map.title ' CDF'],'interpreter','none') @@ -620,7 +620,7 @@ if iload==0 if nfitpost_deciles(jt)) & (y0<=post_deciles(jt+1))); leg{jt}=[int2str(jt) '-dec']; end -[proba] = stab_map_1(x0, indy{1}, indy{end}, [], fname, options_, parnames, estim_params_,0); +[proba] = gsa.stability_mapping_univariate(x0, indy{1}, indy{end}, [], fname, options_, parnames, estim_params_,0); indmcf=find(proba alpha2) under the stable and unacceptable subsets % -% USES qmc_sequence, stab_map_1, stab_map_2 +% USES qmc_sequence, gsa.stability_mapping_univariate, gsa.stability_mapping_bivariate % % Written by Marco Ratto % Joint Research Centre, The European Commission, @@ -147,7 +147,7 @@ if fload==0 %run new MC yys=zeros(length(dr_.ys),Nsam); if opt_gsa.morris == 1 - [lpmat] = Sampling_Function_2(nliv, np+nshock, ntra, ones(np+nshock, 1), zeros(np+nshock,1), []); + [lpmat] = gsa.Sampling_Function_2(nliv, np+nshock, ntra, ones(np+nshock, 1), zeros(np+nshock,1), []); lpmat = lpmat.*(nliv-1)/nliv+1/nliv/2; Nsam=size(lpmat,1); lpmat0 = lpmat(:,1:nshock); @@ -167,7 +167,7 @@ if fload==0 %run new MC end end end - prior_draw_gsa(M_,bayestopt_,options_,estim_params_,1); %initialize + gsa.prior_draw(M_,bayestopt_,options_,estim_params_,1); %initialize if pprior for j=1:nshock if opt_gsa.morris~=1 @@ -184,7 +184,7 @@ if fload==0 %run new MC lpmat(:,j)=lpmat(:,j).*(upper_bound-lower_bound)+lower_bound; end else - xx=prior_draw_gsa(M_,bayestopt_,options_,estim_params_,0,[lpmat0 lpmat]); + xx=gsa.prior_draw(M_,bayestopt_,options_,estim_params_,0,[lpmat0 lpmat]); lpmat0=xx(:,1:nshock); lpmat=xx(:,nshock+1:end); clear xx; @@ -500,7 +500,7 @@ if ~isempty(iunstable) || ~isempty(iwrong) options_mcf.nobeha_title_latex = 'NO unique Stable Saddle-Path'; end options_mcf.title = 'unique solution'; - mcf_analysis(lpmat, istable, itmp, options_mcf, M_, options_, bayestopt_, estim_params_); + gsa.monte_carlo_filtering_analysis(lpmat, istable, itmp, options_mcf, M_, options_, bayestopt_, estim_params_); if ~isempty(iindeterm) itmp = isolve(~ismember(isolve,iindeterm)); @@ -513,7 +513,7 @@ if ~isempty(iunstable) || ~isempty(iwrong) options_mcf.nobeha_title_latex = 'indeterminacy'; end options_mcf.title = 'indeterminacy'; - mcf_analysis(lpmat, itmp, iindeterm, options_mcf, M_, options_, bayestopt_, estim_params_); + gsa.monte_carlo_filtering_analysis(lpmat, itmp, iindeterm, options_mcf, M_, options_, bayestopt_, estim_params_); end if ~isempty(ixun) @@ -527,7 +527,7 @@ if ~isempty(iunstable) || ~isempty(iwrong) options_mcf.nobeha_title_latex = 'explosive solution'; end options_mcf.title = 'instability'; - mcf_analysis(lpmat, itmp, ixun, options_mcf, M_, options_, bayestopt_, estim_params_); + gsa.monte_carlo_filtering_analysis(lpmat, itmp, ixun, options_mcf, M_, options_, bayestopt_, estim_params_); end inorestriction = istable(~ismember(istable,irestriction)); % violation of prior restrictions @@ -543,7 +543,7 @@ if ~isempty(iunstable) || ~isempty(iwrong) options_mcf.nobeha_title_latex = 'inability to find a solution'; end options_mcf.title = 'inability to find a solution'; - mcf_analysis(lpmat, itmp, iwrong, options_mcf, M_, options_, bayestopt_, estim_params_); + gsa.monte_carlo_filtering_analysis(lpmat, itmp, iwrong, options_mcf, M_, options_, bayestopt_, estim_params_); end if ~isempty(irestriction) @@ -576,7 +576,7 @@ if ~isempty(iunstable) || ~isempty(iwrong) options_mcf.nobeha_title_latex = 'NO prior IRF/moment calibration'; end options_mcf.title = 'prior restrictions'; - mcf_analysis([lpmat0 lpmat], irestriction, inorestriction, options_mcf, M_, options_, bayestopt_, estim_params_); + gsa.monte_carlo_filtering_analysis([lpmat0 lpmat], irestriction, inorestriction, options_mcf, M_, options_, bayestopt_, estim_params_); iok = irestriction(1); x0 = [lpmat0(iok,:)'; lpmat(iok,:)']; else diff --git a/matlab/gsa/stab_map_2.m b/matlab/+gsa/stability_mapping_bivariate.m similarity index 96% rename from matlab/gsa/stab_map_2.m rename to matlab/+gsa/stability_mapping_bivariate.m index c6c06da19..4d508c6a4 100644 --- a/matlab/gsa/stab_map_2.m +++ b/matlab/+gsa/stability_mapping_bivariate.m @@ -1,5 +1,5 @@ -function indcorr = stab_map_2(x,alpha2, pvalue_crit, M_,options_,bayestopt_,estim_params_, case_name_plain, case_name_latex, dirname,xparam1,figtitle,fig_caption_latex) -% indcorr = stab_map_2(x,alpha2, pvalue_crit, M_,options_,bayestopt_,estim_params_, fnam, fnam_latex, dirname,xparam1,figtitle,fig_caption_latex) +function indcorr = stability_mapping_bivariate(x,alpha2, pvalue_crit, M_,options_,bayestopt_,estim_params_, case_name_plain, case_name_latex, dirname,xparam1,figtitle,fig_caption_latex) +% indcorr = stability_mapping_bivariate(x,alpha2, pvalue_crit, M_,options_,bayestopt_,estim_params_, fnam, fnam_latex, dirname,xparam1,figtitle,fig_caption_latex) % Inputs: % - x % - alpha2 diff --git a/matlab/gsa/stab_map_1.m b/matlab/+gsa/stability_mapping_univariate.m similarity index 88% rename from matlab/gsa/stab_map_1.m rename to matlab/+gsa/stability_mapping_univariate.m index d6e6d1680..56c9d00c9 100644 --- a/matlab/gsa/stab_map_1.m +++ b/matlab/+gsa/stability_mapping_univariate.m @@ -1,5 +1,5 @@ -function [proba, dproba] = stab_map_1(lpmat, ibehaviour, inonbehaviour, aname, fname_, options_, parnames, estim_params_, iplot, ipar, dirname, pcrit, atitle) -% [proba, dproba] = stab_map_1(lpmat, ibehaviour, inonbehaviour, aname, fname_, options_, parnames, estim_params_, iplot, ipar, dirname, pcrit, atitle) +function [proba, dproba] = stability_mapping_univariate(lpmat, ibehaviour, inonbehaviour, aname, fname_, options_, parnames, estim_params_, iplot, ipar, dirname, pcrit, atitle) +% [proba, dproba] = stability_mapping_univariate(lpmat, ibehaviour, inonbehaviour, aname, fname_, options_, parnames, estim_params_, iplot, ipar, dirname, pcrit, atitle) % Inputs: % - lpmat [double] Monte Carlo matrix % - ibehaviour [integer] index of behavioural runs @@ -18,7 +18,7 @@ function [proba, dproba] = stab_map_1(lpmat, ibehaviour, inonbehaviour, aname, f % % Plots: dotted lines for BEHAVIOURAL % solid lines for NON BEHAVIOURAL -% USES smirnov +% USES gsa.smirnov_test.m % % Written by Marco Ratto % Joint Research Centre, The European Commission, @@ -71,7 +71,7 @@ end proba=NaN(npar,1); dproba=NaN(npar,1); for j=1:npar - [~,P,KSSTAT] = smirnov(lpmat(ibehaviour,j),lpmat(inonbehaviour,j)); + [~,P,KSSTAT] = gsa.smirnov_test(lpmat(ibehaviour,j),lpmat(inonbehaviour,j)); proba(j)=P; dproba(j)=KSSTAT; end @@ -88,12 +88,12 @@ if iplot && ~options_.nograph for j=1+12*(i-1):min(nparplot,12*i) subplot(3,4,j-12*(i-1)) if ~isempty(ibehaviour) - h=cumplot(lpmat(ibehaviour,j)); + h=gsa.cumplot(lpmat(ibehaviour,j)); set(h,'color',[0 0 1], 'linestyle',':','LineWidth',1.5) end hold on if ~isempty(inonbehaviour) - h=cumplot(lpmat(inonbehaviour,j)); + h=gsa.cumplot(lpmat(inonbehaviour,j)); set(h,'color',[0 0 0],'LineWidth',1.5) end title([ftit{j},'. p-value ', num2str(proba(ipar(j)),2)],'interpreter','none') @@ -102,7 +102,7 @@ if iplot && ~options_.nograph dyn_saveas(hh_fig,[dirname,filesep,fname_,'_',aname,'_SA_',int2str(i)],options_.nodisplay,options_.graph_format); if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format))) fidTeX = fopen([dirname,filesep,fname_,'_',aname,'_SA_',int2str(i) '.tex'],'w'); - fprintf(fidTeX,'%% TeX eps-loader file generated by stab_map_1.m (Dynare).\n'); + fprintf(fidTeX,'%% TeX eps-loader file generated by gsa.stability_mapping_univariate.m (Dynare).\n'); fprintf(fidTeX,['%% ' datestr(now,0) '\n\n']); fprintf(fidTeX,'\\begin{figure}[H]\n'); fprintf(fidTeX,'\\centering \n'); @@ -117,7 +117,7 @@ if iplot && ~options_.nograph dyn_saveas(hh_fig,[dirname,filesep,fname_,'_',aname,'_SA'],options_.nodisplay,options_.graph_format); if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format))) fidTeX = fopen([dirname,filesep,fname_,'_',aname,'_SA.tex'],'w'); - fprintf(fidTeX,'%% TeX eps-loader file generated by stab_map_1.m (Dynare).\n'); + fprintf(fidTeX,'%% TeX eps-loader file generated by gsa.stability_mapping_univariate.m (Dynare).\n'); fprintf(fidTeX,['%% ' datestr(now,0) '\n\n']); fprintf(fidTeX,'\\begin{figure}[H]\n'); fprintf(fidTeX,'\\centering \n'); diff --git a/matlab/gsa/stand_.m b/matlab/+gsa/standardize_columns.m similarity index 93% rename from matlab/gsa/stand_.m rename to matlab/+gsa/standardize_columns.m index c95ccf143..9f59f20c6 100644 --- a/matlab/gsa/stand_.m +++ b/matlab/+gsa/standardize_columns.m @@ -1,5 +1,5 @@ -function [y, meany, stdy] = stand_(x) -% [y, meany, stdy] = stand_(x) +function [y, meany, stdy] = standardize_columns(x) +% [y, meany, stdy] = standardize_columns(x) % Standardise a matrix by columns % % [x,my,sy]=stand(y) diff --git a/matlab/gsa/tcrit.m b/matlab/+gsa/tcrit.m similarity index 100% rename from matlab/gsa/tcrit.m rename to matlab/+gsa/tcrit.m diff --git a/matlab/gsa/teff.m b/matlab/+gsa/teff.m similarity index 100% rename from matlab/gsa/teff.m rename to matlab/+gsa/teff.m diff --git a/matlab/gsa/th_moments.m b/matlab/+gsa/th_moments.m similarity index 100% rename from matlab/gsa/th_moments.m rename to matlab/+gsa/th_moments.m diff --git a/matlab/identification_analysis.m b/matlab/+identification/analysis.m similarity index 94% rename from matlab/identification_analysis.m rename to matlab/+identification/analysis.m index 09c9cefe3..e2dba28ac 100644 --- a/matlab/identification_analysis.m +++ b/matlab/+identification/analysis.m @@ -1,5 +1,5 @@ -function [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide_dynamic, derivatives_info, info, error_indicator] = identification_analysis(M_,options_,oo_,bayestopt_,estim_params_,params, indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, init) -% [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide_dynamic, derivatives_info, info, error_indicator] = identification_analysis(M_,options_,oo_,bayestopt_,estim_params_,params, indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, init) +function [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide_dynamic, derivatives_info, info, error_indicator] = analysis(M_,options_,oo_,bayestopt_,estim_params_,params, indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, init) +% [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide_dynamic, derivatives_info, info, error_indicator] = analysis(M_,options_,oo_,bayestopt_,estim_params_,params, indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, init) % ------------------------------------------------------------------------- % This function wraps all identification analysis, i.e. it % (1) wraps functions for the theoretical identification analysis based on moments (Iskrev, 2010), @@ -58,18 +58,18 @@ function [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide % indicator on problems % ------------------------------------------------------------------------- % This function is called by -% * dynare_identification.m +% * identification.run % ------------------------------------------------------------------------- % This function calls % * [M_.fname,'.dynamic'] % * dseries % * dsge_likelihood.m % * dyn_vech -% * ident_bruteforce -% * identification_checks -% * identification_checks_via_subsets +% * identification.bruteforce +% * identification.checks +% * identification.checks_via_subsets % * isoctave -% * get_identification_jacobians (previously getJJ) +% * identification.get_jacobians (previously getJJ) % * matlab_ver_less_than % * prior_bounds % * resol @@ -120,7 +120,7 @@ if ~isempty(estim_params_) M_ = set_all_parameters(params,estim_params_,M_); end -%get options (see dynare_identification.m for description of options) +%get options (see identification.run.m for description of options) nlags = options_ident.ar; advanced = options_ident.advanced; replic = options_ident.replic; @@ -142,7 +142,7 @@ error_indicator.identification_spectrum=0; if info(1) == 0 %no errors in solution % Compute parameter Jacobians for identification analysis - [~, ~, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dSPECTRUM_NO_MEAN, dMINIMAL, derivatives_info] = get_identification_jacobians(estim_params_, M_, options_, options_ident, indpmodel, indpstderr, indpcorr, indvobs, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state); + [~, ~, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dSPECTRUM_NO_MEAN, dMINIMAL, derivatives_info] = identification.get_jacobians(estim_params_, M_, options_, options_ident, indpmodel, indpstderr, indpcorr, indvobs, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state); if isempty(dMINIMAL) % Komunjer and Ng is not computed if (1) minimality conditions are not fullfilled or (2) there are more shocks and measurement errors than observables, so we need to reset options error_indicator.identification_minimal = 1; @@ -206,7 +206,7 @@ if info(1) == 0 %no errors in solution options_ident_local.no_identification_spectrum = 1; %do not recompute dSPECTRUM options_ident_local.ar = nlags; %store new lag number options_.ar = nlags; %store new lag number - [~, ~, ~, ~, ~, ~, MOMENTS, dMOMENTS, ~, ~, ~, ~] = get_identification_jacobians(estim_params_, M_, options_, options_ident_local, indpmodel, indpstderr, indpcorr, indvobs, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state); + [~, ~, ~, ~, ~, ~, MOMENTS, dMOMENTS, ~, ~, ~, ~] = identification.get_jacobians(estim_params_, M_, options_, options_ident_local, indpmodel, indpstderr, indpcorr, indvobs, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state); ind_dMOMENTS = (find(max(abs(dMOMENTS'),[],1) > tol_deriv)); %new index with non-zero rows end @@ -305,7 +305,7 @@ if info(1) == 0 %no errors in solution options_.analytic_derivation = analytic_derivation; %reset option AHess = -AHess; %take negative of hessian if min(eig(AHess))<-tol_rank - error('identification_analysis: Analytic Hessian is not positive semi-definite!') + error('identification.analysis: Analytic Hessian is not positive semi-definite!') end ide_hess.AHess = AHess; %store asymptotic Hessian %normalize asymptotic hessian @@ -313,9 +313,9 @@ if info(1) == 0 %no errors in solution iflag = any((deltaM.*deltaM)==0); %check if all second-order derivatives wrt to a single parameter are nonzero tildaM = AHess./((deltaM)*(deltaM')); %this normalization is for numerical purposes if iflag || rank(AHess)>rank(tildaM) - [ide_hess.cond, ide_hess.rank, ide_hess.ind0, ide_hess.indno, ide_hess.ino, ide_hess.Mco, ide_hess.Pco] = identification_checks(AHess, 0, tol_rank, tol_sv, totparam_nbr); + [ide_hess.cond, ide_hess.rank, ide_hess.ind0, ide_hess.indno, ide_hess.ino, ide_hess.Mco, ide_hess.Pco] = identification.checks(AHess, 0, tol_rank, tol_sv, totparam_nbr); else %use normalized version if possible - [ide_hess.cond, ide_hess.rank, ide_hess.ind0, ide_hess.indno, ide_hess.ino, ide_hess.Mco, ide_hess.Pco] = identification_checks(tildaM, 0, tol_rank, tol_sv, totparam_nbr); + [ide_hess.cond, ide_hess.rank, ide_hess.ind0, ide_hess.indno, ide_hess.ino, ide_hess.Mco, ide_hess.Pco] = identification.checks(tildaM, 0, tol_rank, tol_sv, totparam_nbr); end indok = find(max(ide_hess.indno,[],1)==0); ide_uncert_unnormaliz(indok) = sqrt(diag(inv(AHess(indok,indok))))'; @@ -325,7 +325,7 @@ if info(1) == 0 %no errors in solution diag_chh = sum(si_dREDUCEDFORM(:,ind1)'.*temp1)'; ind1 = ind1(ind1>stderrparam_nbr+corrparam_nbr); cdynamic = si_dDYNAMIC(:,ind1-stderrparam_nbr-corrparam_nbr)*((AHess(ind1,ind1))\si_dDYNAMIC(:,ind1-stderrparam_nbr-corrparam_nbr)'); - flag_score = 1; %this is used for the title in plot_identification.m + flag_score = 1; %this is used for the title in identification.plot.m catch %Asymptotic Hessian via simulation if options_.order > 1 @@ -336,7 +336,7 @@ if info(1) == 0 %no errors in solution options_.periods = periods+100; end replic = max([replic, length(ind_dMOMENTS)*3]); - cmm = simulated_moment_uncertainty(ind_dMOMENTS, periods, replic,options_,M_,oo_); %covariance matrix of moments + cmm = identification.simulated_moment_uncertainty(ind_dMOMENTS, periods, replic,options_,M_,oo_); %covariance matrix of moments sd = sqrt(diag(cmm)); cc = cmm./(sd*sd'); [VV,DD,WW] = eig(cc); @@ -350,9 +350,9 @@ if info(1) == 0 %no errors in solution iflag = any((deltaM.*deltaM)==0); tildaM = MIM./((deltaM)*(deltaM')); if iflag || rank(MIM)>rank(tildaM) - [ide_hess.cond, ide_hess.rank, ide_hess.ind0, ide_hess.indno, ide_hess.ino, ide_hess.Mco, ide_hess.Pco] = identification_checks(MIM, 0, tol_rank, tol_sv, totparam_nbr); + [ide_hess.cond, ide_hess.rank, ide_hess.ind0, ide_hess.indno, ide_hess.ino, ide_hess.Mco, ide_hess.Pco] = identification.checks(MIM, 0, tol_rank, tol_sv, totparam_nbr); else %use normalized version if possible - [ide_hess.cond, ide_hess.rank, ide_hess.ind0, ide_hess.indno, ide_hess.ino, ide_hess.Mco, ide_hess.Pco] = identification_checks(tildaM, 0, tol_rank, tol_sv, totparam_nbr); + [ide_hess.cond, ide_hess.rank, ide_hess.ind0, ide_hess.indno, ide_hess.ino, ide_hess.Mco, ide_hess.Pco] = identification.checks(tildaM, 0, tol_rank, tol_sv, totparam_nbr); end indok = find(max(ide_hess.indno,[],1)==0); ind1 = find(ide_hess.ind0); @@ -363,7 +363,7 @@ if info(1) == 0 %no errors in solution if ~isempty(indok) ide_uncert_unnormaliz(indok) = (sqrt(diag(inv(tildaM(indok,indok))))./deltaM(indok))'; %sqrt(diag(inv(MIM(indok,indok))))'; end - flag_score = 0; %this is used for the title in plot_identification.m + flag_score = 0; %this is used for the title in identification.plot.m end % end of computing sample information matrix for identification strength measure ide_strength_dMOMENTS(indok) = (1./(ide_uncert_unnormaliz(indok)'./abs(params(indok)'))); %this is s_i in Ratto and Iskrev (2011, p.13) @@ -465,11 +465,11 @@ if info(1) == 0 %no errors in solution ide_moments.MOMENTS = MOMENTS; if advanced - % here we do not normalize (i.e. we set norm_dMOMENTS=1) as the OLS in ident_bruteforce is very sensitive to norm_dMOMENTS - [ide_moments.pars, ide_moments.cosndMOMENTS] = ident_bruteforce(M_.dname,M_.fname,dMOMENTS(ind_dMOMENTS,:), max_dim_cova_group, options_.TeX, options_ident.name_tex, options_ident.tittxt, tol_deriv); + % here we do not normalize (i.e. we set norm_dMOMENTS=1) as the OLS in identification.bruteforce is very sensitive to norm_dMOMENTS + [ide_moments.pars, ide_moments.cosndMOMENTS] = identification.bruteforce(M_.dname,M_.fname,dMOMENTS(ind_dMOMENTS,:), max_dim_cova_group, options_.TeX, options_ident.name_tex, options_ident.tittxt, tol_deriv); end - %here we focus on the unnormalized S and V, which is then used in plot_identification.m and for prior_mc > 1 + %here we focus on the unnormalized S and V, which is then used in identification.plot.m and for prior_mc > 1 [~, S, V] = svd(dMOMENTS(ind_dMOMENTS,:),0); if size(S,1) == 1 S = S(1); % edge case that S is not a matrix but a row vector @@ -522,9 +522,9 @@ if info(1) == 0 %no errors in solution %% Perform identification checks, i.e. find out which parameters are involved if checks_via_subsets - % identification_checks_via_subsets is only for debugging + % identification.checks_via_subsets is only for debugging [ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal] = ... - identification_checks_via_subsets(ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, totparam_nbr, modparam_nbr, options_ident, error_indicator); + identification.checks_via_subsets(ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, totparam_nbr, modparam_nbr, options_ident, error_indicator); if ~error_indicator.identification_minimal ide_minimal.minimal_state_space=1; else @@ -532,19 +532,19 @@ if info(1) == 0 %no errors in solution end else [ide_dynamic.cond, ide_dynamic.rank, ide_dynamic.ind0, ide_dynamic.indno, ide_dynamic.ino, ide_dynamic.Mco, ide_dynamic.Pco, ide_dynamic.jweak, ide_dynamic.jweak_pair] = ... - identification_checks(dDYNAMIC(ind_dDYNAMIC,:)./norm_dDYNAMIC, 1, tol_rank, tol_sv, modparam_nbr); + identification.checks(dDYNAMIC(ind_dDYNAMIC,:)./norm_dDYNAMIC, 1, tol_rank, tol_sv, modparam_nbr); if ~options_ident.no_identification_reducedform && ~error_indicator.identification_reducedform [ide_reducedform.cond, ide_reducedform.rank, ide_reducedform.ind0, ide_reducedform.indno, ide_reducedform.ino, ide_reducedform.Mco, ide_reducedform.Pco, ide_reducedform.jweak, ide_reducedform.jweak_pair] = ... - identification_checks(dREDUCEDFORM(ind_dREDUCEDFORM,:)./norm_dREDUCEDFORM, 1, tol_rank, tol_sv, totparam_nbr); + identification.checks(dREDUCEDFORM(ind_dREDUCEDFORM,:)./norm_dREDUCEDFORM, 1, tol_rank, tol_sv, totparam_nbr); end if ~options_ident.no_identification_moments && ~error_indicator.identification_moments [ide_moments.cond, ide_moments.rank, ide_moments.ind0, ide_moments.indno, ide_moments.ino, ide_moments.Mco, ide_moments.Pco, ide_moments.jweak, ide_moments.jweak_pair] = ... - identification_checks(dMOMENTS(ind_dMOMENTS,:)./norm_dMOMENTS, 1, tol_rank, tol_sv, totparam_nbr); + identification.checks(dMOMENTS(ind_dMOMENTS,:)./norm_dMOMENTS, 1, tol_rank, tol_sv, totparam_nbr); end if ~options_ident.no_identification_minimal if ~error_indicator.identification_minimal [ide_minimal.cond, ide_minimal.rank, ide_minimal.ind0, ide_minimal.indno, ide_minimal.ino, ide_minimal.Mco, ide_minimal.Pco, ide_minimal.jweak, ide_minimal.jweak_pair] = ... - identification_checks(dMINIMAL(ind_dMINIMAL,:)./norm_dMINIMAL, 2, tol_rank, tol_sv, totparam_nbr); + identification.checks(dMINIMAL(ind_dMINIMAL,:)./norm_dMINIMAL, 2, tol_rank, tol_sv, totparam_nbr); ide_minimal.minimal_state_space=1; else ide_minimal.minimal_state_space=0; @@ -552,7 +552,7 @@ if info(1) == 0 %no errors in solution end if ~options_ident.no_identification_spectrum && ~error_indicator.identification_spectrum [ide_spectrum.cond, ide_spectrum.rank, ide_spectrum.ind0, ide_spectrum.indno, ide_spectrum.ino, ide_spectrum.Mco, ide_spectrum.Pco, ide_spectrum.jweak, ide_spectrum.jweak_pair] = ... - identification_checks(tilda_dSPECTRUM, 3, tol_rank, tol_sv, totparam_nbr); + identification.checks(tilda_dSPECTRUM, 3, tol_rank, tol_sv, totparam_nbr); end end end diff --git a/matlab/ident_bruteforce.m b/matlab/+identification/bruteforce.m similarity index 98% rename from matlab/ident_bruteforce.m rename to matlab/+identification/bruteforce.m index 75229b4e8..c4be89e40 100644 --- a/matlab/ident_bruteforce.m +++ b/matlab/+identification/bruteforce.m @@ -18,7 +18,7 @@ function [pars, cosnJ] = ident_bruteforce(dname,fname,J, max_dim_cova_group, TeX % cosnJ : cosn of each column with the selected group of columns % ------------------------------------------------------------------------- % This function is called by -% * identification_analysis.m +% * identification.analysis.m % ========================================================================= % Copyright © 2009-2023 Dynare Team % @@ -67,7 +67,7 @@ for ll = 1:max_dim_cova_group cosnJ2=zeros(size(tmp2,1),1); b=[]; for jj = 1:size(tmp2,1) - [cosnJ2(jj,1), b(:,jj)] = cosn([J(:,ii),J(:,tmp2(jj,:))]); + [cosnJ2(jj,1), b(:,jj)] = identification.cosn([J(:,ii),J(:,tmp2(jj,:))]); end cosnJ(ii,ll) = max(cosnJ2(:,1)); if cosnJ(ii,ll)>tol_deriv diff --git a/matlab/identification_checks.m b/matlab/+identification/checks.m similarity index 95% rename from matlab/identification_checks.m rename to matlab/+identification/checks.m index be54d1be1..71c62c012 100644 --- a/matlab/identification_checks.m +++ b/matlab/+identification/checks.m @@ -1,5 +1,5 @@ -function [condX, rankX, ind0, indno, ixno, Mco, Pco, jweak, jweak_pair] = identification_checks(X, test_flag, tol_rank, tol_sv, param_nbr) -% function [condX, rankX, ind0, indno, ixno, Mco, Pco, jweak, jweak_pair] = identification_checks(X, test_flag, tol_rank, tol_sv, param_nbr) +function [condX, rankX, ind0, indno, ixno, Mco, Pco, jweak, jweak_pair] = checks(X, test_flag, tol_rank, tol_sv, param_nbr) +% function [condX, rankX, ind0, indno, ixno, Mco, Pco, jweak, jweak_pair] = checks(X, test_flag, tol_rank, tol_sv, param_nbr) % ------------------------------------------------------------------------- % Checks rank criteria of identification tests and finds out parameter sets % that are not identifiable via the nullspace, pairwise correlation @@ -24,10 +24,10 @@ function [condX, rankX, ind0, indno, ixno, Mco, Pco, jweak, jweak_pair] = identi % * jweak_pair [(vech) matrix] gives 1 if a couple parameters has Pco=1 (with tolerance tol_rank) % ------------------------------------------------------------------------- % This function is called by -% * identification_analysis.m +% * identification.analysis.m % ------------------------------------------------------------------------- % This function calls -% * cosn +% * identification.cosn % * dyn_vech % * vnorm % ========================================================================= @@ -141,7 +141,7 @@ if test_flag == 0 || test_flag == 3 % G is a Gram matrix and hence should be a c else Mco = NaN(param_nbr,1); for ii = 1:size(Xparnonzero,2) - Mco(ind1(ii),:) = cosn([Xparnonzero(:,ii) , Xparnonzero(:,find([1:1:size(Xparnonzero,2)]~=ii)), Xrest]); + Mco(ind1(ii),:) = identification.cosn([Xparnonzero(:,ii) , Xparnonzero(:,find([1:1:size(Xparnonzero,2)]~=ii)), Xrest]); end end @@ -176,7 +176,7 @@ if test_flag ~= 0 for ii = 1:size(Xparnonzero,2) Pco(ind1(ii),ind1(ii)) = 1; for jj = ii+1:size(Xparnonzero,2) - Pco(ind1(ii),ind1(jj)) = cosn([Xparnonzero(:,ii),Xparnonzero(:,jj),Xrest]); + Pco(ind1(ii),ind1(jj)) = identification.cosn([Xparnonzero(:,ii),Xparnonzero(:,jj),Xrest]); Pco(ind1(jj),ind1(ii)) = Pco(ind1(ii),ind1(jj)); end end diff --git a/matlab/identification_checks_via_subsets.m b/matlab/+identification/checks_via_subsets.m similarity index 98% rename from matlab/identification_checks_via_subsets.m rename to matlab/+identification/checks_via_subsets.m index 871b88242..4d75e3736 100644 --- a/matlab/identification_checks_via_subsets.m +++ b/matlab/+identification/checks_via_subsets.m @@ -1,5 +1,5 @@ -function [ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal] = identification_checks_via_subsets(ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, totparam_nbr, modparam_nbr, options_ident,error_indicator) -%[ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal] = identification_checks_via_subsets(ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, totparam_nbr, modparam_nbr, options_ident,error_indicator) +function [ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal] = checks_via_subsets(ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, totparam_nbr, modparam_nbr, options_ident,error_indicator) +%[ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal] = checks_via_subsets(ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, totparam_nbr, modparam_nbr, options_ident,error_indicator) % ------------------------------------------------------------------------- % Finds problematic sets of paramters via checking the necessary rank condition % of the Jacobians for all possible combinations of parameters. The rank is @@ -50,7 +50,7 @@ function [ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal] % * rank: [integer] rank of Jacobian % ------------------------------------------------------------------------- % This function is called by -% * identification_analysis.m +% * identification.analysis.m % ========================================================================= % Copyright © 2019-2021 Dynare Team % @@ -161,7 +161,7 @@ end % initialize for spectrum criteria if ~no_identification_spectrum && ~error_indicator.identification_spectrum - dSPECTRUM = ide_spectrum.tilda_dSPECTRUM; %tilda dSPECTRUM is normalized dSPECTRUM matrix in identification_analysis.m + dSPECTRUM = ide_spectrum.tilda_dSPECTRUM; %tilda dSPECTRUM is normalized dSPECTRUM matrix in identification.analysis.m %alternative normalization %dSPECTRUM = ide_spectrum.dSPECTRUM; %dSPECTRUM(ide_spectrum.ind_dSPECTRUM,:) = dSPECTRUM(ide_spectrum.ind_dSPECTRUM,:)./ide_spectrum.norm_dSPECTRUM; %normalize diff --git a/matlab/cosn.m b/matlab/+identification/cosn.m similarity index 98% rename from matlab/cosn.m rename to matlab/+identification/cosn.m index 7ccd1b5be..a662c245e 100644 --- a/matlab/cosn.m +++ b/matlab/+identification/cosn.m @@ -17,7 +17,7 @@ function [co, b, yhat] = cosn(H) % * y [n by 1] predicted endogenous values given ols estimation % ------------------------------------------------------------------------- % This function is called by -% * identification_checks.m +% * identification.checks.m % * ident_bruteforce.m % ========================================================================= % Copyright © 2008-2019 Dynare Team diff --git a/matlab/disp_identification.m b/matlab/+identification/display.m similarity index 98% rename from matlab/disp_identification.m rename to matlab/+identification/display.m index c72387432..a0726b868 100644 --- a/matlab/disp_identification.m +++ b/matlab/+identification/display.m @@ -1,5 +1,5 @@ -function disp_identification(pdraws, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, name, options_ident) -% disp_identification(pdraws, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, name, options_ident) +function display(pdraws, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, name, options_ident) +% display(pdraws, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, name, options_ident) % ------------------------------------------------------------------------- % This function displays all identification analysis to the command line % ========================================================================= @@ -26,7 +26,7 @@ function disp_identification(pdraws, ide_reducedform, ide_moments, ide_spectrum, % * all output is printed on the command line % ------------------------------------------------------------------------- % This function is called by -% * dynare_identification.m +% * identification.run % ========================================================================= % Copyright © 2010-2021 Dynare Team % @@ -207,7 +207,7 @@ for jide = 1:4 end end - %% display problematic parameters computed by identification_checks_via_subsets + %% display problematic parameters computed by identification.checks_via_subsets elseif checks_via_subsets if ide.rank < size(Jacob,2) no_warning_message_display = 0; diff --git a/matlab/get_identification_jacobians.m b/matlab/+identification/get_jacobians.m similarity index 97% rename from matlab/get_identification_jacobians.m rename to matlab/+identification/get_jacobians.m index c60c0d1af..fc1ba1436 100644 --- a/matlab/get_identification_jacobians.m +++ b/matlab/+identification/get_jacobians.m @@ -1,5 +1,5 @@ -function [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dSPECTRUM_NO_MEAN, dMINIMAL, derivatives_info] = get_identification_jacobians(estim_params, M_, options_, options_ident, indpmodel, indpstderr, indpcorr, indvobs, dr, endo_steady_state, exo_steady_state, exo_det_steady_state) -% [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dSPECTRUM_NO_MEAN, dMINIMAL, derivatives_info] = get_identification_jacobians(estim_params, M_, options_, options_ident, indpmodel, indpstderr, indpcorr, indvobs, dr, endo_steady_state, exo_steady_state, exo_det_steady_state) +function [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dSPECTRUM_NO_MEAN, dMINIMAL, derivatives_info] = get_jacobians(estim_params, M_, options_, options_ident, indpmodel, indpstderr, indpcorr, indvobs, dr, endo_steady_state, exo_steady_state, exo_det_steady_state) +% [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dSPECTRUM_NO_MEAN, dMINIMAL, derivatives_info] = get_jacobians(estim_params, M_, options_, options_ident, indpmodel, indpstderr, indpcorr, indvobs, dr, endo_steady_state, exo_steady_state, exo_det_steady_state) % previously getJJ.m in Dynare 4.5 % Sets up the Jacobians needed for identification analysis % ========================================================================= @@ -84,7 +84,7 @@ function [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dM % % ------------------------------------------------------------------------- % This function is called by -% * identification_analysis.m +% * identification.analysis.m % ------------------------------------------------------------------------- % This function calls % * commutation @@ -94,7 +94,7 @@ function [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dM % * fjaco % * get_perturbation_params_derivs (previously getH) % * get_all_parameters -% * identification_numerical_objective (previously thet2tau) +% * identification.numerical_objective (previously thet2tau) % * pruned_state_space_system % * vec % ========================================================================= @@ -258,7 +258,7 @@ if ~no_identification_moments if kronflag == -1 %numerical derivative of autocovariogram - dMOMENTS = fjaco(str2func('identification_numerical_objective'), xparam1, 1, estim_params, M_, options_, indpmodel, indpstderr, indvobs, useautocorr, nlags, grid_nbr, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); %[outputflag=1] + dMOMENTS = fjaco(str2func('identification.numerical_objective'), xparam1, 1, estim_params, M_, options_, indpmodel, indpstderr, indvobs, useautocorr, nlags, grid_nbr, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); %[outputflag=1] dMOMENTS = [dMEAN; dMOMENTS]; %add Jacobian of steady state of VAROBS variables else dMOMENTS = zeros(obs_nbr + obs_nbr*(obs_nbr+1)/2 + nlags*obs_nbr^2 , totparam_nbr); @@ -315,7 +315,7 @@ if ~no_identification_spectrum IA = eye(size(pruned.A,1)); if kronflag == -1 %numerical derivative of spectral density - dOmega_tmp = fjaco(str2func('identification_numerical_objective'), xparam1, 2, estim_params, M_, options_, indpmodel, indpstderr, indvobs, useautocorr, nlags, grid_nbr, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); %[outputflag=2] + dOmega_tmp = fjaco(str2func('identification.numerical_objective'), xparam1, 2, estim_params, M_, options_, indpmodel, indpstderr, indvobs, useautocorr, nlags, grid_nbr, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); %[outputflag=2] kk = 0; for ig = 1:length(freqs) kk = kk+1; diff --git a/matlab/identification_numerical_objective.m b/matlab/+identification/numerical_objective.m similarity index 97% rename from matlab/identification_numerical_objective.m rename to matlab/+identification/numerical_objective.m index 548e68784..ec0ffc793 100644 --- a/matlab/identification_numerical_objective.m +++ b/matlab/+identification/numerical_objective.m @@ -22,7 +22,7 @@ function out = identification_numerical_objective(params, outputflag, estim_para % 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. @@ -32,7 +32,7 @@ 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'] diff --git a/matlab/plot_identification.m b/matlab/+identification/plot.m similarity index 95% rename from matlab/plot_identification.m rename to matlab/+identification/plot.m index 035c9a325..d20531f75 100644 --- a/matlab/plot_identification.m +++ b/matlab/+identification/plot.m @@ -1,5 +1,5 @@ -function plot_identification(M_, params, idemoments, idehess, idemodel, idelre, advanced, tittxt, name, IdentifDirectoryName, fname, options_, estim_params_, bayestopt_, tit_TeX, name_tex) -% plot_identification(M_, params,idemoments,idehess,idemodel, idelre, advanced, tittxt, name, IdentifDirectoryName, fname, options_, estim_params_, bayestopt_, tit_TeX, name_tex) +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 @@ -156,7 +156,7 @@ if SampleSize == 1 end if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format))) fidTeX = fopen([IdentifDirectoryName '/' fname '_ident_strength_' tittxt1,'.tex'],'w'); - fprintf(fidTeX,'%% TeX eps-loader file generated by plot_identification.m (Dynare).\n'); + 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'); @@ -203,7 +203,7 @@ if SampleSize == 1 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 '/' fname '_sensitivity_' tittxt1,'.tex'],'w'); - fprintf(fidTeX,'%% TeX eps-loader file generated by plot_identification.m (Dynare).\n'); + 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'); @@ -262,7 +262,7 @@ if SampleSize == 1 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 '/' fname '_ident_collinearity_' tittxt1 '_' int2str(j),'.tex'],'w'); - fprintf(fidTeX,'%% TeX eps-loader file generated by plot_identification.m (Dynare).\n'); + 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'); @@ -329,7 +329,7 @@ if SampleSize == 1 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 '/' fname '_ident_pattern_' tittxt1 '_1','.tex'],'w'); - fprintf(fidTeX,'%% TeX eps-loader file generated by plot_identification.m (Dynare).\n'); + 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'); @@ -344,7 +344,7 @@ if SampleSize == 1 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 '/' fname '_ident_pattern_' tittxt1 '_2.tex'],'w'); - fprintf(fidTeX,'%% TeX eps-loader file generated by plot_identification.m (Dynare).\n'); + 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'); @@ -392,7 +392,7 @@ else 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 '/' fname '_MC_sensitivity.tex'],'w'); - fprintf(fidTeX,'%% TeX eps-loader file generated by plot_identification.m (Dynare).\n'); + 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'); @@ -450,17 +450,17 @@ else options_mcf.title = 'MC Highest Condition Number LRE Model'; ncut=floor(SampleSize/10*9); [~,is]=sort(idelre.cond); - mcf_analysis(params, is(1:ncut), is(ncut+1:end), options_mcf, M_, options_, bayestopt_, estim_params_); + 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'; [~,is]=sort(idemodel.cond); - mcf_analysis(params, is(1:ncut), is(ncut+1:end), options_mcf, M_, options_, bayestopt_, estim_params_); + 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'; [~,is]=sort(idemoments.cond); - mcf_analysis(params, is(1:ncut), is(ncut+1:end), options_mcf, M_, options_, bayestopt_, estim_params_); + 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']); @@ -514,7 +514,7 @@ else 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 '/' fname '_MC_ident_pattern_1.tex'],'w'); - fprintf(fidTeX,'%% TeX eps-loader file generated by plot_identification.m (Dynare).\n'); + 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'); @@ -529,7 +529,7 @@ else 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 '/' fname '_MC_ident_pattern_2.tex'],'w'); - fprintf(fidTeX,'%% TeX eps-loader file generated by plot_identification.m (Dynare).\n'); + 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'); diff --git a/matlab/dynare_identification.m b/matlab/+identification/run.m similarity index 94% rename from matlab/dynare_identification.m rename to matlab/+identification/run.m index 1d9d23dd8..e716fa1ea 100644 --- a/matlab/dynare_identification.m +++ b/matlab/+identification/run.m @@ -1,5 +1,5 @@ -function [pdraws, STO_REDUCEDFORM, STO_MOMENTS, STO_DYNAMIC, STO_si_dDYNAMIC, STO_si_dREDUCEDFORM, STO_si_dMOMENTS, STO_dSPECTRUM, STO_dMINIMAL] = dynare_identification(M_,oo_,options_,bayestopt_,estim_params_,options_ident, pdraws0) -% [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 @@ -32,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 @@ -95,7 +95,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) @@ -140,7 +140,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. @@ -181,7 +181,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; @@ -308,12 +308,12 @@ options_.options_ident = []; 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 @@ -384,7 +384,7 @@ else % no estimated_params block, choose all model parameters and all stderr par 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; @@ -402,13 +402,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 @@ -471,7 +471,7 @@ if iload <=0 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(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 + 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_); @@ -488,7 +488,7 @@ if iload <=0 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(M_,options_,oo_,bayestopt_,estim_params_,params, indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1); + identification.analysis(M_,options_,oo_,bayestopt_,estim_params_,params, indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1); end end if info(1) @@ -513,10 +513,10 @@ 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(M_,params, ide_moments_point, ide_hess_point, ide_reducedform_point, ide_dynamic_point, options_ident.advanced, parameters, name, ... + 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 @@ -529,7 +529,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 @@ -543,7 +543,7 @@ if iload <=0 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(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 + 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 @@ -801,25 +801,25 @@ if iload <=0 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); 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; 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; end if ~options_MC.no_identification_spectrum && ~isempty(STO_dSPECTRUM) - % note that this is not the same dSPECTRUMnorm as computed in identification_analysis + % note that this is not the same dSPECTRUMnorm as computed in identification.analysis dSPECTRUMnorm(iter,:) = vnorm(STO_dSPECTRUM(:,:,irun)); %not yet used end if ~options_MC.no_identification_minimal && ~isempty(STO_dMINIMAL) - % note that this is not the same dMINIMALnorm as computed in identification_analysis + % note that this is not the same dMINIMALnorm as computed in identification.analysis dMINIMALnorm(iter,:) = vnorm(STO_dMINIMAL(:,:,irun)); %not yet used end end @@ -847,7 +847,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 = []; @@ -876,10 +876,10 @@ 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(M_,ide_hess_point.params, ide_moments_point, ide_hess_point, ide_reducedform_point, ide_dynamic_point, options_ident.advanced, parameters, name, ... + 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 @@ -890,11 +890,11 @@ 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(M_, pdraws, IDE_MOMENTS, ide_hess_point, IDE_REDUCEDFORM, IDE_DYNAMIC, options_ident.advanced, 'MC sample ', name, ... + 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 @@ -912,15 +912,15 @@ if SampleSize > 1 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(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 + 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(M_, pdraws(jmax,:), ide_moments_max, ide_hess_max, ide_reducedform_max, ide_dynamic_max, 1, tittxt, name, ... + 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 @@ -931,15 +931,15 @@ if SampleSize > 1 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, ~, ~, 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 + 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(M_, pdraws(jmin,:),ide_moments_min,ide_hess_min,ide_reducedform_min,ide_dynamic_min,1,tittxt,name,... + 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 @@ -954,14 +954,14 @@ 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(M_,options_,oo_,bayestopt_,estim_params_,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(M_, pdraws(jcrit(j),:), ide_moments_(j), ide_hess_(j), ide_reducedform_(j), ide_dynamic_(j), 1, tittxt, name, ... + 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 diff --git a/matlab/simulated_moment_uncertainty.m b/matlab/+identification/simulated_moment_uncertainty.m similarity index 100% rename from matlab/simulated_moment_uncertainty.m rename to matlab/+identification/simulated_moment_uncertainty.m diff --git a/matlab/commutation.m b/matlab/commutation.m index f0c8c6aa5..1fffe85e8 100644 --- a/matlab/commutation.m +++ b/matlab/commutation.m @@ -14,7 +14,7 @@ function k = commutation(n, m, sparseflag) % ------------------------------------------------------------------------- % This function is called by % * get_perturbation_params_derivs.m (previously getH.m) -% * get_identification_jacobians.m (previously getJJ.m) +% * identification.get_jacobians.m (previously getJJ.m) % * pruned_state_space_system.m % ------------------------------------------------------------------------- % This function calls diff --git a/matlab/duplication.m b/matlab/duplication.m index 9afecf520..c69a6719f 100644 --- a/matlab/duplication.m +++ b/matlab/duplication.m @@ -11,7 +11,7 @@ function [Dp,DpMPinv] = duplication(p) % DpMPinv: Moore-Penroze inverse of Dp % ------------------------------------------------------------------------- % This function is called by -% * get_identification_jacobians.m (previously getJJ.m) +% * identification.get_jacobians.m (previously getJJ.m) % ========================================================================= % Copyright © 1997 Tom Minka % Copyright © 2019 Dynare Team diff --git a/matlab/fjaco.m b/matlab/fjaco.m index 3d41787b1..b020ed269 100644 --- a/matlab/fjaco.m +++ b/matlab/fjaco.m @@ -30,7 +30,7 @@ function fjac = fjaco(f,x,varargin) ff=feval(f,x,varargin{:}); tol = eps.^(1/3); %some default value -if strcmp(func2str(f),'get_perturbation_params_derivs_numerical_objective') || strcmp(func2str(f),'identification_numerical_objective') +if strcmp(func2str(f),'get_perturbation_params_derivs_numerical_objective') || strcmp(func2str(f),'identification.numerical_objective') tol= varargin{4}.dynatol.x; end h = tol.*max(abs(x),1); @@ -40,12 +40,12 @@ fjac = NaN(length(ff),length(x)); for j=1:length(x) xx = x; xx(j) = xh1(j); f1=feval(f,xx,varargin{:}); - if isempty(f1) && (strcmp(func2str(f),'get_perturbation_params_derivs_numerical_objective') || strcmp(func2str(f),'identification_numerical_objective') ) + if isempty(f1) && (strcmp(func2str(f),'get_perturbation_params_derivs_numerical_objective') || strcmp(func2str(f),'identification.numerical_objective') ) [~,info]=feval(f,xx,varargin{:}); disp_info_error_identification_perturbation(info,j); end xx(j) = xh0(j); f0=feval(f,xx,varargin{:}); - if isempty(f0) && (strcmp(func2str(f),'get_perturbation_params_derivs_numerical_objective') || strcmp(func2str(f),'identification_numerical_objective') ) + if isempty(f0) && (strcmp(func2str(f),'get_perturbation_params_derivs_numerical_objective') || strcmp(func2str(f),'identification.numerical_objective') ) [~,info]=feval(f,xx,varargin{:}); disp_info_error_identification_perturbation(info,j) end diff --git a/matlab/get_minimal_state_representation.m b/matlab/get_minimal_state_representation.m index 277717f3b..deb6d43bb 100644 --- a/matlab/get_minimal_state_representation.m +++ b/matlab/get_minimal_state_representation.m @@ -53,7 +53,7 @@ function [CheckCO,minns,minSYS] = get_minimal_state_representation(SYS, derivs_f % Jacobian (wrt to all parameters) of measurement matrix minD % ------------------------------------------------------------------------- % This function is called by -% * get_identification_jacobians.m (previously getJJ.m) +% * identification.get_jacobians.m (previously getJJ.m) % ------------------------------------------------------------------------- % This function calls % * check_minimality (embedded) diff --git a/matlab/get_perturbation_params_derivs.m b/matlab/get_perturbation_params_derivs.m index fab20ef03..6721b3a42 100644 --- a/matlab/get_perturbation_params_derivs.m +++ b/matlab/get_perturbation_params_derivs.m @@ -88,7 +88,7 @@ function DERIVS = get_perturbation_params_derivs(M_, options_, estim_params_, dr % ------------------------------------------------------------------------- % This function is called by % * dsge_likelihood.m -% * get_identification_jacobians.m +% * identification.get_jacobians.m % ------------------------------------------------------------------------- % This function calls % * [fname,'.dynamic'] diff --git a/matlab/list_of_functions_to_be_cleared.m b/matlab/list_of_functions_to_be_cleared.m index 78a18484b..200478bb2 100644 --- a/matlab/list_of_functions_to_be_cleared.m +++ b/matlab/list_of_functions_to_be_cleared.m @@ -1,2 +1,2 @@ -list_of_functions = {'discretionary_policy_1', 'dsge_var_likelihood', 'dyn_first_order_solver', 'dyn_waitbar', 'ep_residuals', 'evaluate_likelihood', 'prior_draw_gsa', 'identification_analysis', 'computeDLIK', 'univariate_computeDLIK', 'metropolis_draw', 'flag_implicit_skip_nan', 'mr_hessian', 'masterParallel', 'auxiliary_initialization', 'auxiliary_particle_filter', 'conditional_filter_proposal', 'conditional_particle_filter', 'gaussian_filter', 'gaussian_filter_bank', 'gaussian_mixture_filter', 'gaussian_mixture_filter_bank', 'Kalman_filter', 'online_auxiliary_filter', 'pruned_state_space_system', 'sequential_importance_particle_filter', 'solve_model_for_online_filter', 'prior_draw', 'priordens',... +list_of_functions = {'discretionary_policy_1', 'dsge_var_likelihood', 'dyn_first_order_solver', 'dyn_waitbar', 'ep_residuals', 'evaluate_likelihood', '+gsa/prior_draw.m', '+identification/analysis.m', 'computeDLIK', 'univariate_computeDLIK', 'metropolis_draw', 'flag_implicit_skip_nan', 'mr_hessian', 'masterParallel', 'auxiliary_initialization', 'auxiliary_particle_filter', 'conditional_filter_proposal', 'conditional_particle_filter', 'gaussian_filter', 'gaussian_filter_bank', 'gaussian_mixture_filter', 'gaussian_mixture_filter_bank', 'Kalman_filter', 'online_auxiliary_filter', 'pruned_state_space_system', 'sequential_importance_particle_filter', 'solve_model_for_online_filter', 'prior_draw', 'priordens',... '+occbin/solver.m','+occbin/mkdatap_anticipated_dyn.m','+occbin/mkdatap_anticipated_2constraints_dyn.m','+occbin/match_function.m','+occbin/solve_one_constraint.m','+occbin/solve_two_constraint.m','+occbin/plot/shock_decomposition.m'}; diff --git a/matlab/pruned_state_space_system.m b/matlab/pruned_state_space_system.m index 3f1f51a14..40d974ca9 100644 --- a/matlab/pruned_state_space_system.m +++ b/matlab/pruned_state_space_system.m @@ -80,8 +80,8 @@ function pruned_state_space = pruned_state_space_system(M_, options_, dr, indy, % parameter Jacobian of E_y % ------------------------------------------------------------------------- % This function is called by -% * get_identification_jacobians.m -% * identification_numerical_objective.m +% * identification.get_jacobians.m +% * identification.numerical_objective.m % ------------------------------------------------------------------------- % This function calls % * allVL1.m diff --git a/matlab/set_all_parameters.m b/matlab/set_all_parameters.m index c97f98aa2..ed090a176 100644 --- a/matlab/set_all_parameters.m +++ b/matlab/set_all_parameters.m @@ -26,7 +26,7 @@ function M_ = set_all_parameters(xparam1,estim_params_,M_) %! @sp 2 %! @strong{This function is called by:} %! @sp 1 -%! @ref{DsgeSmoother}, @ref{dynare_estimation_1}, @ref{@@gsa/filt_mc_}, @ref{identification_analysis}, @ref{PosteriorFilterSmootherAndForecast}, @ref{prior_posterior_statistics_core}, @ref{prior_sampler} +%! @ref{DsgeSmoother}, @ref{dynare_estimation_1}, @ref{@@gsa.monte_carlo_filtering}, @ref{identification.analysis}, @ref{PosteriorFilterSmootherAndForecast}, @ref{prior_posterior_statistics_core}, @ref{prior_sampler} %! @sp 2 %! @strong{This function calls:} %! @sp 2 diff --git a/preprocessor b/preprocessor index 3dadac8f1..8d0e8cca5 160000 --- a/preprocessor +++ b/preprocessor @@ -1 +1 @@ -Subproject commit 3dadac8f191dfa1dd660442d5bc4526c4c218149 +Subproject commit 8d0e8cca5cb78b9dde0ecc867ffb0c64d06dd338 From c3268c02795eb75d235ee425dd4a7186ed4501d2 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Fri, 8 Dec 2023 10:06:00 +0100 Subject: [PATCH 2/3] Move various functions from main matlab folder to subfolders --- license.txt | 8 +- matlab/+identification/analysis.m | 6 +- matlab/+identification/checks.m | 2 +- matlab/{ => +identification}/fjaco.m | 6 +- matlab/+identification/get_jacobians.m | 14 +-- .../get_minimal_state_representation.m | 0 .../get_perturbation_params_derivs.m | 28 +++--- ...bation_params_derivs_numerical_objective.m | 2 +- matlab/{ => +identification}/hessian_sparse.m | 0 matlab/+identification/numerical_objective.m | 6 +- matlab/+identification/run.m | 10 +- matlab/{ => +identification}/unfold_g3.m | 0 matlab/{ => +identification}/unfold_g4.m | 0 matlab/{ => +identification}/vnorm.m | 0 matlab/+mom/objective_function.m | 6 +- matlab/+osr/objective.m | 4 +- matlab/+pruned_SS/Q6_plication.m | 2 +- matlab/{ => +pruned_SS}/commutation.m | 0 matlab/{ => +pruned_SS}/duplication.m | 0 .../pruned_state_space_system.m | 12 +-- matlab/+pruned_SS/quadruplication.m | 2 +- matlab/{ => +pruned_SS}/uperm.m | 96 +++++++++---------- matlab/{ => backward}/simul_static_model.m | 0 matlab/cellofchar2mfile.m | 63 ------------ matlab/check_model.m | 2 +- matlab/clear_persistent_variables.m | 50 +++++++++- .../{ => cli}/print_moments_implied_prior.m | 0 matlab/compute_model_moments.m | 32 ------- matlab/convergence_diagnostics/mcmc_ifac.m | 24 +++++ matlab/dyn_autocorr.m | 40 -------- matlab/dyn_diag_vech.m | 31 ------ matlab/dynare_config.m | 6 ++ matlab/dynare_gradient.m | 66 ------------- matlab/{ => ep}/cartesian_product_of_sets.m | 0 .../{ => ep}/cubature_with_gaussian_weight.m | 0 .../gauss_hermite_weights_and_nodes.m | 0 .../gauss_legendre_weights_and_nodes.m | 0 matlab/estimation/dsge_likelihood.m | 4 +- .../dsge_simulated_theoretical_covariance.m | 2 +- matlab/evaluate_dynamic_model.m | 29 ------ matlab/{ => kalman}/DsgeSmoother.m | 0 matlab/{ => kalman}/evaluate_smoother.m | 0 .../missing_DiffuseKalmanSmootherH1_Z.m | 0 .../missing_DiffuseKalmanSmootherH3_Z.m | 0 .../save_display_classical_smoother_results.m | 0 matlab/{ => kalman}/store_smoother_results.m | 0 matlab/{ => latex}/collect_latex_files.m | 0 matlab/{ => latex}/dyn_latex_table.m | 0 matlab/{ => latex}/isbayes.m | 0 matlab/{ => latex}/write_latex_definitions.m | 0 .../{ => latex}/write_latex_parameter_table.m | 0 matlab/{ => latex}/write_latex_prior_table.m | 0 matlab/{ => lmmcp}/mcp_func.m | 0 matlab/{ => matrix_solver}/fastgensylv.m | 0 matlab/{ => matrix_solver}/gensylv_fp.m | 0 matlab/{ => matrix_solver}/lyapunov_solver.m | 0 matlab/{ => matrix_solver}/lyapunov_symm.m | 0 .../quadratic_matrix_equation_solver.m | 0 matlab/{ => matrix_solver}/sylvester3.m | 0 matlab/{ => matrix_solver}/sylvester3a.m | 0 .../{ => moments}/UnivariateSpectralDensity.m | 0 matlab/{ => moments}/add_filter_subtitle.m | 0 .../{ => moments}/compute_moments_varendo.m | 0 .../conditional_variance_decomposition.m | 0 ...al_variance_decomposition_ME_mc_analysis.m | 0 ...ional_variance_decomposition_mc_analysis.m | 0 .../{ => moments}/correlation_mc_analysis.m | 0 matlab/{ => moments}/covariance_mc_analysis.m | 0 matlab/{ => moments}/disp_moments.m | 0 .../disp_th_moments_pruned_state_space.m | 2 +- matlab/my_subplot.m | 56 ----------- .../reduced_rank_cholesky.m | 0 .../{ => optimal_policy}/dyn_ramsey_static.m | 0 .../evaluate_planner_objective.m | 0 .../get_optimal_policy_discount_factor.m | 0 .../{ => optimal_policy}/mult_elimination.m | 0 matlab/{ => optimization}/csolve.m | 0 matlab/{ => optimization}/dynare_solve.m | 0 matlab/{ => optimization}/lnsrch1.m | 0 matlab/{ => optimization}/options2cell.m | 0 matlab/{ => optimization}/solve1.m | 0 .../{ => optimization}/solve_one_boundary.m | 0 .../step_length_correction.m | 0 matlab/{ => optimization}/trust_region.m | 0 .../basic_plan.m | 0 .../evaluate_max_dynamic_residual.m | 0 .../flip_plan.m | 0 .../init_plan.m | 0 matlab/save_results.m | 38 -------- .../WriteShockDecomp2Excel.m | 0 .../annualized_shock_decomposition.m | 0 .../epilogue_shock_decomposition.m | 0 .../{ => shock_decomposition}/expand_group.m | 0 .../{ => shock_decomposition}/graph_decomp.m | 0 .../graph_decomp_detail.m | 0 .../initial_condition_decomposition.m | 0 .../plot_shock_decomposition.m | 0 .../realtime_shock_decomposition.m | 0 ..._initial_condition_decomposition_options.m | 0 ...default_plot_shock_decomposition_options.m | 0 .../shock_decomposition.m | 0 .../squeeze_shock_decomposition.m | 0 .../AIM_first_order_solver.m | 0 .../compute_decision_rules.m | 0 .../convertAimCodeToInfo.m | 0 matlab/{ => stochastic_solver}/disp_dr.m | 0 .../dyn_first_order_solver.m | 0 .../dyn_second_order_solver.m | 0 .../{ => stochastic_solver}/dynare_resolve.m | 0 .../getIrfShocksIndx.m | 0 matlab/{ => stochastic_solver}/irf.m | 0 matlab/{ => stochastic_solver}/k_order_pert.m | 0 .../kalman_transition_matrix.m | 0 matlab/{ => stochastic_solver}/resol.m | 0 .../select_qz_criterium_value.m | 0 .../{ => stochastic_solver}/set_state_space.m | 0 matlab/{ => stochastic_solver}/simult.m | 0 matlab/{ => stochastic_solver}/simult_.m | 0 matlab/{ => stochastic_solver}/simultxdet.m | 0 matlab/{ => stochastic_solver}/stoch_simul.m | 0 .../stochastic_solvers.m | 0 .../BrockMirman_PertParamsDerivs.mod | 4 +- .../burnside_3_order_PertParamsDerivs.mod | 2 +- .../as2007_minimal.mod | 2 +- .../minimal_state_space_system/sw_minimal.mod | 2 +- .../AS_pruned_state_space_red_shock.mod | 2 +- .../AnSchorfheide_pruned_state_space.mod | 2 +- .../pruning/fs2000_pruning.mod | 4 +- 128 files changed, 195 insertions(+), 472 deletions(-) rename matlab/{ => +identification}/fjaco.m (83%) rename matlab/{ => +identification}/get_minimal_state_representation.m (100%) rename matlab/{ => +identification}/get_perturbation_params_derivs.m (97%) rename matlab/{ => +identification}/get_perturbation_params_derivs_numerical_objective.m (98%) rename matlab/{ => +identification}/hessian_sparse.m (100%) rename matlab/{ => +identification}/unfold_g3.m (100%) rename matlab/{ => +identification}/unfold_g4.m (100%) rename matlab/{ => +identification}/vnorm.m (100%) rename matlab/{ => +pruned_SS}/commutation.m (100%) rename matlab/{ => +pruned_SS}/duplication.m (100%) rename matlab/{ => +pruned_SS}/pruned_state_space_system.m (99%) rename matlab/{ => +pruned_SS}/uperm.m (96%) rename matlab/{ => backward}/simul_static_model.m (100%) delete mode 100644 matlab/cellofchar2mfile.m rename matlab/{ => cli}/print_moments_implied_prior.m (100%) delete mode 100644 matlab/compute_model_moments.m delete mode 100644 matlab/dyn_autocorr.m delete mode 100644 matlab/dyn_diag_vech.m delete mode 100644 matlab/dynare_gradient.m rename matlab/{ => ep}/cartesian_product_of_sets.m (100%) rename matlab/{ => ep}/cubature_with_gaussian_weight.m (100%) rename matlab/{ => ep}/gauss_hermite_weights_and_nodes.m (100%) rename matlab/{ => ep}/gauss_legendre_weights_and_nodes.m (100%) delete mode 100644 matlab/evaluate_dynamic_model.m rename matlab/{ => kalman}/DsgeSmoother.m (100%) rename matlab/{ => kalman}/evaluate_smoother.m (100%) rename matlab/{ => kalman}/missing_DiffuseKalmanSmootherH1_Z.m (100%) rename matlab/{ => kalman}/missing_DiffuseKalmanSmootherH3_Z.m (100%) rename matlab/{ => kalman}/save_display_classical_smoother_results.m (100%) rename matlab/{ => kalman}/store_smoother_results.m (100%) rename matlab/{ => latex}/collect_latex_files.m (100%) rename matlab/{ => latex}/dyn_latex_table.m (100%) rename matlab/{ => latex}/isbayes.m (100%) rename matlab/{ => latex}/write_latex_definitions.m (100%) rename matlab/{ => latex}/write_latex_parameter_table.m (100%) rename matlab/{ => latex}/write_latex_prior_table.m (100%) rename matlab/{ => lmmcp}/mcp_func.m (100%) rename matlab/{ => matrix_solver}/fastgensylv.m (100%) rename matlab/{ => matrix_solver}/gensylv_fp.m (100%) rename matlab/{ => matrix_solver}/lyapunov_solver.m (100%) rename matlab/{ => matrix_solver}/lyapunov_symm.m (100%) rename matlab/{ => matrix_solver}/quadratic_matrix_equation_solver.m (100%) rename matlab/{ => matrix_solver}/sylvester3.m (100%) rename matlab/{ => matrix_solver}/sylvester3a.m (100%) rename matlab/{ => moments}/UnivariateSpectralDensity.m (100%) rename matlab/{ => moments}/add_filter_subtitle.m (100%) rename matlab/{ => moments}/compute_moments_varendo.m (100%) rename matlab/{ => moments}/conditional_variance_decomposition.m (100%) rename matlab/{ => moments}/conditional_variance_decomposition_ME_mc_analysis.m (100%) rename matlab/{ => moments}/conditional_variance_decomposition_mc_analysis.m (100%) rename matlab/{ => moments}/correlation_mc_analysis.m (100%) rename matlab/{ => moments}/covariance_mc_analysis.m (100%) rename matlab/{ => moments}/disp_moments.m (100%) delete mode 100644 matlab/my_subplot.m rename matlab/{ => nonlinear-filters}/reduced_rank_cholesky.m (100%) rename matlab/{ => optimal_policy}/dyn_ramsey_static.m (100%) rename matlab/{ => optimal_policy}/evaluate_planner_objective.m (100%) rename matlab/{ => optimal_policy}/get_optimal_policy_discount_factor.m (100%) rename matlab/{ => optimal_policy}/mult_elimination.m (100%) rename matlab/{ => optimization}/csolve.m (100%) rename matlab/{ => optimization}/dynare_solve.m (100%) rename matlab/{ => optimization}/lnsrch1.m (100%) rename matlab/{ => optimization}/options2cell.m (100%) rename matlab/{ => optimization}/solve1.m (100%) rename matlab/{ => optimization}/solve_one_boundary.m (100%) rename matlab/{ => optimization}/step_length_correction.m (100%) rename matlab/{ => optimization}/trust_region.m (100%) rename matlab/{ => perfect-foresight-models}/basic_plan.m (100%) rename matlab/{ => perfect-foresight-models}/evaluate_max_dynamic_residual.m (100%) rename matlab/{ => perfect-foresight-models}/flip_plan.m (100%) rename matlab/{ => perfect-foresight-models}/init_plan.m (100%) delete mode 100644 matlab/save_results.m rename matlab/{ => shock_decomposition}/WriteShockDecomp2Excel.m (100%) rename matlab/{ => shock_decomposition}/annualized_shock_decomposition.m (100%) rename matlab/{ => shock_decomposition}/epilogue_shock_decomposition.m (100%) rename matlab/{ => shock_decomposition}/expand_group.m (100%) rename matlab/{ => shock_decomposition}/graph_decomp.m (100%) rename matlab/{ => shock_decomposition}/graph_decomp_detail.m (100%) rename matlab/{ => shock_decomposition}/initial_condition_decomposition.m (100%) rename matlab/{ => shock_decomposition}/plot_shock_decomposition.m (100%) rename matlab/{ => shock_decomposition}/realtime_shock_decomposition.m (100%) rename matlab/{ => shock_decomposition}/set_default_initial_condition_decomposition_options.m (100%) rename matlab/{ => shock_decomposition}/set_default_plot_shock_decomposition_options.m (100%) rename matlab/{ => shock_decomposition}/shock_decomposition.m (100%) rename matlab/{ => shock_decomposition}/squeeze_shock_decomposition.m (100%) rename matlab/{ => stochastic_solver}/AIM_first_order_solver.m (100%) rename matlab/{ => stochastic_solver}/compute_decision_rules.m (100%) rename matlab/{ => stochastic_solver}/convertAimCodeToInfo.m (100%) rename matlab/{ => stochastic_solver}/disp_dr.m (100%) rename matlab/{ => stochastic_solver}/dyn_first_order_solver.m (100%) rename matlab/{ => stochastic_solver}/dyn_second_order_solver.m (100%) rename matlab/{ => stochastic_solver}/dynare_resolve.m (100%) rename matlab/{ => stochastic_solver}/getIrfShocksIndx.m (100%) rename matlab/{ => stochastic_solver}/irf.m (100%) rename matlab/{ => stochastic_solver}/k_order_pert.m (100%) rename matlab/{ => stochastic_solver}/kalman_transition_matrix.m (100%) rename matlab/{ => stochastic_solver}/resol.m (100%) rename matlab/{ => stochastic_solver}/select_qz_criterium_value.m (100%) rename matlab/{ => stochastic_solver}/set_state_space.m (100%) rename matlab/{ => stochastic_solver}/simult.m (100%) rename matlab/{ => stochastic_solver}/simult_.m (100%) rename matlab/{ => stochastic_solver}/simultxdet.m (100%) rename matlab/{ => stochastic_solver}/stoch_simul.m (100%) rename matlab/{ => stochastic_solver}/stochastic_solvers.m (100%) diff --git a/license.txt b/license.txt index c1bee92b0..73fb4bc4d 100644 --- a/license.txt +++ b/license.txt @@ -113,7 +113,7 @@ Copyright: 1995 E.G.Tsionas 2015-2017 Dynare Team License: GPL-3+ -Files: matlab/endogenous_prior.m +Files: matlab/estimation/endogenous_prior.m Copyright: 2011 Lawrence J. Christiano, Mathias Trabandt and Karl Walentin 2013-2017 Dynare Team License: GPL-3+ @@ -128,7 +128,7 @@ Copyright: 2016 Benjamin Born and Johannes Pfeifer 2016-2017 Dynare Team License: GPL-3+ -Files: matlab/commutation.m matlab/duplication.m +Files: matlab/+pruned_SS/commutation.m matlab/+pruned_SS/duplication.m Copyright: 1997 Tom Minka 2019-2020 Dynare Team License: GPL-3+ @@ -141,7 +141,7 @@ Comment: The original author gave authorization to change the license from BSD-2-clause to GPL-3+ and redistribute it under GPL-3+ with Dynare. -Files: matlab/uperm.m +Files: matlab/+pruned_SS/uperm.m Copyright: 2014 Bruno Luong 2020 Dynare Team License: GPL-3+ @@ -149,7 +149,7 @@ Comment: The original author gave authorization to change the license from BSD-2-clause to GPL-3+ and redistribute it under GPL-3+ with Dynare. -Files: matlab/prodmom.m matlab/bivmom.m +Files: matlab/+pruned_SS/prodmom.m matlab/+pruned_SS/bivmom.m Copyright: 2008-2015 Raymond Kan 2019-2020 Dynare Team License: GPL-3+ diff --git a/matlab/+identification/analysis.m b/matlab/+identification/analysis.m index e2dba28ac..c1a387e2a 100644 --- a/matlab/+identification/analysis.m +++ b/matlab/+identification/analysis.m @@ -375,7 +375,7 @@ if info(1) == 0 %no errors in solution if size(quant,1)==1 si_dMOMENTSnorm = abs(quant).*normaliz_prior_std; else - si_dMOMENTSnorm = vnorm(quant).*normaliz_prior_std; + si_dMOMENTSnorm = identification.vnorm(quant).*normaliz_prior_std; end iy = find(diag_chh); ind_dREDUCEDFORM = ind_dREDUCEDFORM(iy); @@ -385,7 +385,7 @@ if info(1) == 0 %no errors in solution if size(quant,1)==1 si_dREDUCEDFORMnorm = abs(quant).*normaliz_prior_std; else - si_dREDUCEDFORMnorm = vnorm(quant).*normaliz_prior_std; + si_dREDUCEDFORMnorm = identification.vnorm(quant).*normaliz_prior_std; end else si_dREDUCEDFORMnorm = []; @@ -399,7 +399,7 @@ if info(1) == 0 %no errors in solution if size(quant,1)==1 si_dDYNAMICnorm = abs(quant).*normaliz_prior_std(stderrparam_nbr+corrparam_nbr+1:end); else - si_dDYNAMICnorm = vnorm(quant).*normaliz_prior_std(stderrparam_nbr+corrparam_nbr+1:end); + si_dDYNAMICnorm = identification.vnorm(quant).*normaliz_prior_std(stderrparam_nbr+corrparam_nbr+1:end); end else si_dDYNAMICnorm=[]; diff --git a/matlab/+identification/checks.m b/matlab/+identification/checks.m index 71c62c012..03895b6ad 100644 --- a/matlab/+identification/checks.m +++ b/matlab/+identification/checks.m @@ -75,7 +75,7 @@ end % find non-zero columns at machine precision if size(Xpar,1) > 1 - ind1 = find(vnorm(Xpar) >= eps); + ind1 = find(identification.vnorm(Xpar) >= eps); else ind1 = find(abs(Xpar) >= eps); % if only one parameter end diff --git a/matlab/fjaco.m b/matlab/+identification/fjaco.m similarity index 83% rename from matlab/fjaco.m rename to matlab/+identification/fjaco.m index b020ed269..b1818932c 100644 --- a/matlab/fjaco.m +++ b/matlab/+identification/fjaco.m @@ -30,7 +30,7 @@ function fjac = fjaco(f,x,varargin) ff=feval(f,x,varargin{:}); tol = eps.^(1/3); %some default value -if strcmp(func2str(f),'get_perturbation_params_derivs_numerical_objective') || strcmp(func2str(f),'identification.numerical_objective') +if strcmp(func2str(f),'identification.get_perturbation_params_derivs_numerical_objective') || strcmp(func2str(f),'identification.numerical_objective') tol= varargin{4}.dynatol.x; end h = tol.*max(abs(x),1); @@ -40,12 +40,12 @@ fjac = NaN(length(ff),length(x)); for j=1:length(x) xx = x; xx(j) = xh1(j); f1=feval(f,xx,varargin{:}); - if isempty(f1) && (strcmp(func2str(f),'get_perturbation_params_derivs_numerical_objective') || strcmp(func2str(f),'identification.numerical_objective') ) + if isempty(f1) && (strcmp(func2str(f),'identification.get_perturbation_params_derivs_numerical_objective') || strcmp(func2str(f),'identification.numerical_objective') ) [~,info]=feval(f,xx,varargin{:}); disp_info_error_identification_perturbation(info,j); end xx(j) = xh0(j); f0=feval(f,xx,varargin{:}); - if isempty(f0) && (strcmp(func2str(f),'get_perturbation_params_derivs_numerical_objective') || strcmp(func2str(f),'identification.numerical_objective') ) + if isempty(f0) && (strcmp(func2str(f),'identification.get_perturbation_params_derivs_numerical_objective') || strcmp(func2str(f),'identification.numerical_objective') ) [~,info]=feval(f,xx,varargin{:}); disp_info_error_identification_perturbation(info,j) end diff --git a/matlab/+identification/get_jacobians.m b/matlab/+identification/get_jacobians.m index fc1ba1436..f1bdac9ef 100644 --- a/matlab/+identification/get_jacobians.m +++ b/matlab/+identification/get_jacobians.m @@ -153,7 +153,7 @@ obs_nbr = length(indvobs); d2flag = 0; % do not compute second parameter derivatives % Get Jacobians (wrt selected params) of steady state, dynamic model derivatives and perturbation solution matrices for all endogenous variables -dr.derivs = get_perturbation_params_derivs(M_, options_, estim_params, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, indpmodel, indpstderr, indpcorr, d2flag); +dr.derivs = identification.get_perturbation_params_derivs(M_, options_, estim_params, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, indpmodel, indpstderr, indpcorr, d2flag); [I,~] = find(lead_lag_incidence'); %I is used to select nonzero columns of the Jacobian of endogenous variables in dynamic model files yy0 = dr.ys(I); %steady state of dynamic (endogenous and auxiliary variables) in lead_lag_incidence order @@ -230,7 +230,7 @@ elseif order == 3 end % Get (pruned) state space representation: -pruned = pruned_state_space_system(M_, options_, dr, indvobs, nlags, useautocorr, 1); +pruned = pruned_SS.pruned_state_space_system(M_, options_, dr, indvobs, nlags, useautocorr, 1); MEAN = pruned.E_y; dMEAN = pruned.dE_y; %storage for Jacobians used in dsge_likelihood.m for analytical Gradient and Hession of likelihood (only at order=1) @@ -258,7 +258,7 @@ if ~no_identification_moments if kronflag == -1 %numerical derivative of autocovariogram - dMOMENTS = fjaco(str2func('identification.numerical_objective'), xparam1, 1, estim_params, M_, options_, indpmodel, indpstderr, indvobs, useautocorr, nlags, grid_nbr, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); %[outputflag=1] + dMOMENTS = identification.fjaco(str2func('identification.numerical_objective'), xparam1, 1, estim_params, M_, options_, indpmodel, indpstderr, indvobs, useautocorr, nlags, grid_nbr, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); %[outputflag=1] dMOMENTS = [dMEAN; dMOMENTS]; %add Jacobian of steady state of VAROBS variables else dMOMENTS = zeros(obs_nbr + obs_nbr*(obs_nbr+1)/2 + nlags*obs_nbr^2 , totparam_nbr); @@ -315,7 +315,7 @@ if ~no_identification_spectrum IA = eye(size(pruned.A,1)); if kronflag == -1 %numerical derivative of spectral density - dOmega_tmp = fjaco(str2func('identification.numerical_objective'), xparam1, 2, estim_params, M_, options_, indpmodel, indpstderr, indvobs, useautocorr, nlags, grid_nbr, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); %[outputflag=2] + dOmega_tmp = identification.fjaco(str2func('identification.numerical_objective'), xparam1, 2, estim_params, M_, options_, indpmodel, indpstderr, indvobs, useautocorr, nlags, grid_nbr, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); %[outputflag=2] kk = 0; for ig = 1:length(freqs) kk = kk+1; @@ -333,7 +333,7 @@ if ~no_identification_spectrum dC = reshape(pruned.dC,size(pruned.dC,1)*size(pruned.dC,2),size(pruned.dC,3)); dD = reshape(pruned.dD,size(pruned.dD,1)*size(pruned.dD,2),size(pruned.dD,3)); dVarinov = reshape(pruned.dVarinov,size(pruned.dVarinov,1)*size(pruned.dVarinov,2),size(pruned.dVarinov,3)); - K_obs_exo = commutation(obs_nbr,size(pruned.Varinov,1)); + K_obs_exo = pruned_SS.commutation(obs_nbr,size(pruned.Varinov,1)); for ig=1:length(freqs) z = tneg(ig); zIminusA = (z*IA - pruned.A); @@ -400,7 +400,7 @@ if ~no_identification_minimal SYS.dC = dr.derivs.dghx(pruned.indy,:,:); SYS.D = dr.ghu(pruned.indy,:); SYS.dD = dr.derivs.dghu(pruned.indy,:,:); - [CheckCO,minnx,SYS] = get_minimal_state_representation(SYS,1); + [CheckCO,minnx,SYS] = identification.get_minimal_state_representation(SYS,1); if CheckCO == 0 warning_KomunjerNg = 'WARNING: Komunjer and Ng (2011) failed:\n'; @@ -423,7 +423,7 @@ if ~no_identification_minimal dvechSig = dvechSig(indvechSig,:); Inx = eye(minnx); Inu = eye(exo_nbr); - [~,Enu] = duplication(exo_nbr); + [~,Enu] = pruned_SS.duplication(exo_nbr); KomunjerNg_DL = [dminA; dminB; dminC; dminD; dvechSig]; KomunjerNg_DT = [kron(transpose(minA),Inx) - kron(Inx,minA); kron(transpose(minB),Inx); diff --git a/matlab/get_minimal_state_representation.m b/matlab/+identification/get_minimal_state_representation.m similarity index 100% rename from matlab/get_minimal_state_representation.m rename to matlab/+identification/get_minimal_state_representation.m diff --git a/matlab/get_perturbation_params_derivs.m b/matlab/+identification/get_perturbation_params_derivs.m similarity index 97% rename from matlab/get_perturbation_params_derivs.m rename to matlab/+identification/get_perturbation_params_derivs.m index 6721b3a42..9f6543998 100644 --- a/matlab/get_perturbation_params_derivs.m +++ b/matlab/+identification/get_perturbation_params_derivs.m @@ -191,7 +191,7 @@ if order > 1 && analytic_derivation_mode == 1 analytic_derivation_mode = 0; fprintf('As order > 1, reset ''analytic_derivation_mode'' to 0\n'); end -numerical_objective_fname = str2func('get_perturbation_params_derivs_numerical_objective'); +numerical_objective_fname = str2func('identification.get_perturbation_params_derivs_numerical_objective'); idx_states = nstatic+(1:nspred); %index for state variables, in DR order modparam_nbr = length(indpmodel); %number of selected model parameters stderrparam_nbr = length(indpstderr); %number of selected stderr parameters @@ -295,7 +295,7 @@ if analytic_derivation_mode == -1 % - perturbation solution matrices: dghx, dghu, dghxx, dghxu, dghuu, dghs2, dghxxx, dghxxu, dghxuu, dghuuu, dghxss, dghuss %Parameter Jacobian of covariance matrix and solution matrices (wrt selected stderr, corr and model paramters) - dSig_gh = fjaco(numerical_objective_fname, xparam1, 'perturbation_solution', estim_params_, M_, options_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); + dSig_gh = identification.fjaco(numerical_objective_fname, xparam1, 'perturbation_solution', estim_params_, M_, options_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); ind_Sigma_e = (1:exo_nbr^2); ind_ghx = ind_Sigma_e(end) + (1:endo_nbr*nspred); ind_ghu = ind_ghx(end) + (1:endo_nbr*exo_nbr); @@ -348,7 +348,7 @@ if analytic_derivation_mode == -1 end %Parameter Jacobian of dynamic model derivatives (wrt selected model parameters only) - dYss_g = fjaco(numerical_objective_fname, modparam1, 'dynamic_model', estim_params_model, M_, options_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); + dYss_g = identification.fjaco(numerical_objective_fname, modparam1, 'dynamic_model', estim_params_model, M_, options_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); ind_Yss = 1:endo_nbr; if options_.discretionary_policy || options_.ramsey_policy ind_g1 = ind_Yss(end) + (1:M_.eq_nbr*yy0ex0_nbr); @@ -374,7 +374,7 @@ if analytic_derivation_mode == -1 % Hessian (wrt paramters) of steady state and first-order solution matrices ghx and Om % note that hessian_sparse.m (contrary to hessian.m) does not take symmetry into account, but focuses already on unique values options_.order = 1; %make sure only first order - d2Yss_KalmanA_Om = hessian_sparse(numerical_objective_fname, xparam1, gstep, 'Kalman_Transition', estim_params_, M_, options_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); + d2Yss_KalmanA_Om = identification.hessian_sparse(numerical_objective_fname, xparam1, gstep, 'Kalman_Transition', estim_params_, M_, options_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); options_.order = order; %make sure to set back ind_KalmanA = ind_Yss(end) + (1:endo_nbr^2); DERIVS.d2KalmanA = d2Yss_KalmanA_Om(ind_KalmanA, indp2tottot2); %only unique elements @@ -394,7 +394,7 @@ if analytic_derivation_mode == -2 % The parameter derivatives of perturbation solution matrices are computed analytically below (analytic_derivation_mode=0) if order == 3 [~, g1, g2, g3] = feval([fname,'.dynamic'], ys(I), exo_steady_state', params, ys, 1); - g3 = unfold_g3(g3, yy0ex0_nbr); + g3 = identification.unfold_g3(g3, yy0ex0_nbr); elseif order == 2 [~, g1, g2] = feval([fname,'.dynamic'], ys(I), exo_steady_state', params, ys, 1); elseif order == 1 @@ -405,7 +405,7 @@ if analytic_derivation_mode == -2 % computation of d2Yss and d2g1 % note that hessian_sparse does not take symmetry into account, i.e. compare hessian_sparse.m to hessian.m, but focuses already on unique values, which are duplicated below options_.order = 1; %d2flag requires only first order - d2Yss_g1 = hessian_sparse(numerical_objective_fname, modparam1, gstep, 'dynamic_model', estim_params_model, M_, options_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); % d2flag requires only first-order + d2Yss_g1 = identification.hessian_sparse(numerical_objective_fname, modparam1, gstep, 'dynamic_model', estim_params_model, M_, options_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); % d2flag requires only first-order options_.order = order; %make sure to set back the order d2Yss = reshape(full(d2Yss_g1(1:endo_nbr,:)), [endo_nbr modparam_nbr modparam_nbr]); %put into tensor notation for j=1:endo_nbr @@ -431,7 +431,7 @@ if analytic_derivation_mode == -2 end %Parameter Jacobian of dynamic model derivatives (wrt selected model parameters only) - dYss_g = fjaco(numerical_objective_fname, modparam1, 'dynamic_model', estim_params_model, M_, options_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); + dYss_g = identification.fjaco(numerical_objective_fname, modparam1, 'dynamic_model', estim_params_model, M_, options_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); ind_Yss = 1:endo_nbr; ind_g1 = ind_Yss(end) + (1:endo_nbr*yy0ex0_nbr); dYss = dYss_g(ind_Yss,:); %in tensor notation, wrt selected model parameters only @@ -460,7 +460,7 @@ elseif (analytic_derivation_mode == 0 || analytic_derivation_mode == 1) [~, ~, g2_static] = feval([fname,'.static'], ys, exo_steady_state', params); %g2_static is [endo_nbr by endo_nbr^2] second derivative (wrt all endogenous variables) of static model equations f, i.e. d(df/dys)/dys, in declaration order if order < 3 [~, g1, g2, g3] = feval([fname,'.dynamic'], ys(I), exo_steady_state', params, ys, 1); %note that g3 does not contain symmetric elements - g3 = unfold_g3(g3, yy0ex0_nbr); %add symmetric elements to g3 + g3 = identification.unfold_g3(g3, yy0ex0_nbr); %add symmetric elements to g3 else T = NaN(sum(dynamic_tmp_nbr(1:5))); T = feval([fname, '.dynamic_g4_tt'], T, ys(I), exo_steady_state', params, ys, 1); @@ -468,8 +468,8 @@ elseif (analytic_derivation_mode == 0 || analytic_derivation_mode == 1) g2 = feval([fname, '.dynamic_g2'], T, ys(I), exo_steady_state', params, ys, 1, false); %g2 is [endo_nbr by yy0ex0_nbr^2] second derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order g3 = feval([fname, '.dynamic_g3'], T, ys(I), exo_steady_state', params, ys, 1, false); %note that g3 does not contain symmetric elements g4 = feval([fname, '.dynamic_g4'], T, ys(I), exo_steady_state', params, ys, 1, false); %note that g4 does not contain symmetric elements - g3 = unfold_g3(g3, yy0ex0_nbr); %add symmetric elements to g3, %g3 is [endo_nbr by yy0ex0_nbr^3] third-derivative (wrt all dynamic variables) of dynamic model equations, i.e. (d(df/dyy0ex0)/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order - g4 = unfold_g4(g4, yy0ex0_nbr); %add symmetric elements to g4, %g4 is [endo_nbr by yy0ex0_nbr^4] fourth-derivative (wrt all dynamic variables) of dynamic model equations, i.e. ((d(df/dyy0ex0)/dyy0ex0)/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order + g3 = identification.unfold_g3(g3, yy0ex0_nbr); %add symmetric elements to g3, %g3 is [endo_nbr by yy0ex0_nbr^3] third-derivative (wrt all dynamic variables) of dynamic model equations, i.e. (d(df/dyy0ex0)/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order + g4 = identification.unfold_g4(g4, yy0ex0_nbr); %add symmetric elements to g4, %g4 is [endo_nbr by yy0ex0_nbr^4] fourth-derivative (wrt all dynamic variables) of dynamic model equations, i.e. ((d(df/dyy0ex0)/dyy0ex0)/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order end %g1 is [endo_nbr by yy0ex0_nbr first derivative (wrt all dynamic variables) of dynamic model equations, i.e. df/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order %g2 is [endo_nbr by yy0ex0_nbr^2] second derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order @@ -556,7 +556,7 @@ elseif (analytic_derivation_mode == 0 || analytic_derivation_mode == 1) error('For analytical parameter derivatives ''dynamic_params_derivs.m'' file is needed, this can be created by putting identification(order=%d) into your mod file.',order) end [~, g1, g2, g3] = feval([fname,'.dynamic'], ys(I), exo_steady_state', params, ys, 1); %note that g3 does not contain symmetric elements - g3 = unfold_g3(g3, yy0ex0_nbr); %add symmetric elements to g3 + g3 = identification.unfold_g3(g3, yy0ex0_nbr); %add symmetric elements to g3 %g1 is [endo_nbr by yy0ex0_nbr first derivative (wrt all dynamic variables) of dynamic model equations, i.e. df/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order %g2 is [endo_nbr by yy0ex0_nbr^2] second derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order %g3 is [endo_nbr by yy0ex0_nbr^3] third-derivative (wrt all dynamic variables) of dynamic model equations, i.e. (d(df/dyy0ex0)/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order @@ -575,8 +575,8 @@ elseif (analytic_derivation_mode == 0 || analytic_derivation_mode == 1) g2 = feval([fname, '.dynamic_g2'], T, ys(I), exo_steady_state', params, ys, 1, false); %g2 is [endo_nbr by yy0ex0_nbr^2] second derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order g3 = feval([fname, '.dynamic_g3'], T, ys(I), exo_steady_state', params, ys, 1, false); %note that g3 does not contain symmetric elements g4 = feval([fname, '.dynamic_g4'], T, ys(I), exo_steady_state', params, ys, 1, false); %note that g4 does not contain symmetric elements - g3 = unfold_g3(g3, yy0ex0_nbr); %add symmetric elements to g3, %g3 is [endo_nbr by yy0ex0_nbr^3] third-derivative (wrt all dynamic variables) of dynamic model equations, i.e. (d(df/dyy0ex0)/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order - g4 = unfold_g4(g4, yy0ex0_nbr); %add symmetric elements to g4, %g4 is [endo_nbr by yy0ex0_nbr^4] fourth-derivative (wrt all dynamic variables) of dynamic model equations, i.e. ((d(df/dyy0ex0)/dyy0ex0)/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order + g3 = identification.unfold_g3(g3, yy0ex0_nbr); %add symmetric elements to g3, %g3 is [endo_nbr by yy0ex0_nbr^3] third-derivative (wrt all dynamic variables) of dynamic model equations, i.e. (d(df/dyy0ex0)/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order + g4 = identification.unfold_g4(g4, yy0ex0_nbr); %add symmetric elements to g4, %g4 is [endo_nbr by yy0ex0_nbr^4] fourth-derivative (wrt all dynamic variables) of dynamic model equations, i.e. ((d(df/dyy0ex0)/dyy0ex0)/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order end end % Parameter Jacobian of steady state in different orderings, note dys is in declaration order @@ -801,7 +801,7 @@ if analytic_derivation_mode == 1 dghu = [zeros(endo_nbr*exo_nbr, stderrparam_nbr+corrparam_nbr) dghu]; % Compute dOm = dvec(ghu*Sigma_e*ghu') from expressions 34 in Iskrev (2010) Appendix A - dOm = kron(I_endo,ghu*Sigma_e)*(commutation(endo_nbr, exo_nbr)*dghu)... + dOm = kron(I_endo,ghu*Sigma_e)*(pruned_SS.commutation(endo_nbr, exo_nbr)*dghu)... + kron(ghu,ghu)*reshape(dSigma_e, exo_nbr^2, totparam_nbr) + kron(ghu*Sigma_e,I_endo)*dghu; % Put into tensor notation diff --git a/matlab/get_perturbation_params_derivs_numerical_objective.m b/matlab/+identification/get_perturbation_params_derivs_numerical_objective.m similarity index 98% rename from matlab/get_perturbation_params_derivs_numerical_objective.m rename to matlab/+identification/get_perturbation_params_derivs_numerical_objective.m index efb997309..4fc7f47f5 100644 --- a/matlab/get_perturbation_params_derivs_numerical_objective.m +++ b/matlab/+identification/get_perturbation_params_derivs_numerical_objective.m @@ -95,7 +95,7 @@ if strcmp(outputflag,'dynamic_model') out = [Yss; g1(:); g2(:)]; elseif options_.order == 3 [~, g1, g2, g3] = feval([M_.fname,'.dynamic'], ys(I), exo_steady_state', M_.params, ys, 1); - g3 = unfold_g3(g3, length(ys(I))+M_.exo_nbr); + g3 = identification.unfold_g3(g3, length(ys(I))+M_.exo_nbr); out = [Yss; g1(:); g2(:); g3(:)]; end end diff --git a/matlab/hessian_sparse.m b/matlab/+identification/hessian_sparse.m similarity index 100% rename from matlab/hessian_sparse.m rename to matlab/+identification/hessian_sparse.m diff --git a/matlab/+identification/numerical_objective.m b/matlab/+identification/numerical_objective.m index ec0ffc793..84c1695b7 100644 --- a/matlab/+identification/numerical_objective.m +++ b/matlab/+identification/numerical_objective.m @@ -1,5 +1,5 @@ -function out = identification_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 = identification_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) +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 @@ -80,7 +80,7 @@ end %% compute Kalman transition matrices and steady state with updated parameters [dr,info,M_.params] = compute_decision_rules(M_,options_,dr, steady_state, exo_steady_state, exo_det_steady_state); options_ = rmfield(options_,'options_ident'); -pruned = pruned_state_space_system(M_, options_, dr, indvar, nlags, useautocorr, 0); +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 diff --git a/matlab/+identification/run.m b/matlab/+identification/run.m index e716fa1ea..fac0c5d22 100644 --- a/matlab/+identification/run.m +++ b/matlab/+identification/run.m @@ -803,24 +803,24 @@ if iload <=0 iter=iter+1; % 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 % 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 % 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 + 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 + dMINIMALnorm(iter,:) = identification.vnorm(STO_dMINIMAL(:,:,irun)); %not yet used end end end diff --git a/matlab/unfold_g3.m b/matlab/+identification/unfold_g3.m similarity index 100% rename from matlab/unfold_g3.m rename to matlab/+identification/unfold_g3.m diff --git a/matlab/unfold_g4.m b/matlab/+identification/unfold_g4.m similarity index 100% rename from matlab/unfold_g4.m rename to matlab/+identification/unfold_g4.m diff --git a/matlab/vnorm.m b/matlab/+identification/vnorm.m similarity index 100% rename from matlab/vnorm.m rename to matlab/+identification/vnorm.m diff --git a/matlab/+mom/objective_function.m b/matlab/+mom/objective_function.m index 42c459856..2e7ed196b 100644 --- a/matlab/+mom/objective_function.m +++ b/matlab/+mom/objective_function.m @@ -150,11 +150,11 @@ if strcmp(options_mom_.mom.mom_method,'GMM') 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_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state, indpmodel, indpstderr, indpcorr, 0); %analytic derivatives of perturbation matrices + oo_.dr.derivs = identification.get_perturbation_params_derivs(M_, options_mom_, estim_params_, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state, 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); + pruned_state_space = pruned_SS.pruned_state_space_system(M_, options_mom_, oo_.dr, oo_.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_, oo_.dr, oo_.mom.obs_var, options_mom_.ar, 0, 0); end oo_.mom.model_moments = NaN(options_mom_.mom.mom_nbr,1); for jm = 1:size(M_.matched_moments,1) diff --git a/matlab/+osr/objective.m b/matlab/+osr/objective.m index 1bb3148ad..5c3b38059 100644 --- a/matlab/+osr/objective.m +++ b/matlab/+osr/objective.m @@ -64,9 +64,9 @@ if ~options_.analytic_derivation loss = full(weights(:)'*vx(:)); else totparam_nbr=length(i_params); - oo_.dr.derivs = get_perturbation_params_derivs(M_, options_, [], oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state, i_params, [], [], 0); %analytic derivatives of perturbation matrices + oo_.dr.derivs = identification.get_perturbation_params_derivs(M_, options_, [], oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state, i_params, [], [], 0); %analytic derivatives of perturbation matrices - pruned_state_space = pruned_state_space_system(M_, options_, oo_.dr, i_var, 0, 0, 1); + pruned_state_space = pruned_SS.pruned_state_space_system(M_, options_, oo_.dr, i_var, 0, 0, 1); vx = pruned_state_space.Var_y + pruned_state_space.E_y*pruned_state_space.E_y'; dE_yy = pruned_state_space.dVar_y; for jp=1:length(i_params) diff --git a/matlab/+pruned_SS/Q6_plication.m b/matlab/+pruned_SS/Q6_plication.m index 68a2bf816..ef432d1c2 100644 --- a/matlab/+pruned_SS/Q6_plication.m +++ b/matlab/+pruned_SS/Q6_plication.m @@ -49,7 +49,7 @@ for i1=1:p for i4=i3:p for i5=i4:p for i6=i5:p - idx = uperm([i6 i5 i4 i3 i2 i1]); + idx = pruned_SS.uperm([i6 i5 i4 i3 i2 i1]); for r = 1:size(idx,1) ii1 = idx(r,1); ii2= idx(r,2); ii3=idx(r,3); ii4=idx(r,4); ii5=idx(r,5); ii6=idx(r,6); n = ii1 + (ii2-1)*p + (ii3-1)*p^2 + (ii4-1)*p^3 + (ii5-1)*p^4 + (ii6-1)*p^5; diff --git a/matlab/commutation.m b/matlab/+pruned_SS/commutation.m similarity index 100% rename from matlab/commutation.m rename to matlab/+pruned_SS/commutation.m diff --git a/matlab/duplication.m b/matlab/+pruned_SS/duplication.m similarity index 100% rename from matlab/duplication.m rename to matlab/+pruned_SS/duplication.m diff --git a/matlab/pruned_state_space_system.m b/matlab/+pruned_SS/pruned_state_space_system.m similarity index 99% rename from matlab/pruned_state_space_system.m rename to matlab/+pruned_SS/pruned_state_space_system.m index 40d974ca9..26d5eeb3d 100644 --- a/matlab/pruned_state_space_system.m +++ b/matlab/+pruned_SS/pruned_state_space_system.m @@ -418,7 +418,7 @@ if order > 1 hx_hu = kron(hx,hu); hu_hu = kron(hu,hu); I_xx = eye(x_nbr^2); - K_x_x = commutation(x_nbr,x_nbr,1); + K_x_x = pruned_SS.commutation(x_nbr,x_nbr,1); invIx_hx = (eye(x_nbr)-hx)\eye(x_nbr); %Compute unique fourth order product moments of u, i.e. unique(E[kron(kron(kron(u,u),u),u)],'stable') @@ -596,9 +596,9 @@ if order > 1 if order > 2 % Some common and useful objects for order > 2 if isempty(K_u_xx) - K_u_xx = commutation(u_nbr,x_nbr^2,1); - K_u_ux = commutation(u_nbr,u_nbr*x_nbr,1); - K_xx_x = commutation(x_nbr^2,x_nbr); + K_u_xx = pruned_SS.commutation(u_nbr,x_nbr^2,1); + K_u_ux = pruned_SS.commutation(u_nbr,u_nbr*x_nbr,1); + K_xx_x = pruned_SS.commutation(x_nbr^2,x_nbr); end hx_hss2 = kron(hx,1/2*hss); hu_hss2 = kron(hu,1/2*hss); @@ -627,7 +627,7 @@ if order > 1 E_xf_xfxs = Var_z(id_z3_xf_xf, id_z2_xs ) + E_xfxf(:)*E_xs'; %this is E[kron(xf,xf)*xs'] E_xf_xfxf_xf = Var_z(id_z3_xf_xf, id_z3_xf_xf) + E_xfxf(:)*E_xfxf(:)'; %this is E[kron(xf,xf)*kron(xf,xf)'] E_xrdxf = reshape(invIxx_hx_hx*vec(... - hxx*reshape( commutation(x_nbr^2,x_nbr,1)*E_xf_xfxs(:), x_nbr^2,x_nbr)*hx'... + hxx*reshape( pruned_SS.commutation(x_nbr^2,x_nbr,1)*E_xf_xfxs(:), x_nbr^2,x_nbr)*hx'... + hxu*kron(E_xs,E_uu)*hu'... + 1/6*hxxx*reshape(E_xf_xfxf_xf,x_nbr^3,x_nbr)*hx'... + 1/6*huuu*reshape(QPu*E_u_u_u_u,u_nbr^3,u_nbr)*hu'... @@ -655,7 +655,7 @@ if order > 1 dE_xf_xfxf_xf(:,:,jp2) = dVar_z(id_z3_xf_xf , id_z3_xf_xf , jp2) + vec(dE_xfxf(:,:,jp2))*E_xfxf(:)' + E_xfxf(:)*vec(dE_xfxf(:,:,jp2))'; dE_xrdxf(:,:,jp2) = reshape(invIxx_hx_hx*vec(... dhx(:,:,jp2)*E_xrdxf*hx' + hx*E_xrdxf*dhx(:,:,jp2)'... - + dhxx(:,:,jp2)*reshape( commutation(x_nbr^2,x_nbr,1)*E_xf_xfxs(:), x_nbr^2,x_nbr)*hx' + hxx*reshape( commutation(x_nbr^2,x_nbr,1)*vec(dE_xf_xfxs(:,:,jp2)), x_nbr^2,x_nbr)*hx' + hxx*reshape( commutation(x_nbr^2,x_nbr,1)*E_xf_xfxs(:), x_nbr^2,x_nbr)*dhx(:,:,jp2)'... + + dhxx(:,:,jp2)*reshape( pruned_SS.commutation(x_nbr^2,x_nbr,1)*E_xf_xfxs(:), x_nbr^2,x_nbr)*hx' + hxx*reshape( pruned_SS.commutation(x_nbr^2,x_nbr,1)*vec(dE_xf_xfxs(:,:,jp2)), x_nbr^2,x_nbr)*hx' + hxx*reshape( pruned_SS.commutation(x_nbr^2,x_nbr,1)*E_xf_xfxs(:), x_nbr^2,x_nbr)*dhx(:,:,jp2)'... + dhxu(:,:,jp2)*kron(E_xs,E_uu)*hu' + hxu*kron(dE_xs(:,jp2),E_uu)*hu' + hxu*kron(E_xs,dE_uu(:,:,jp2))*hu' + hxu*kron(E_xs,E_uu)*dhu(:,:,jp2)'... + 1/6*dhxxx(:,:,jp2)*reshape(E_xf_xfxf_xf,x_nbr^3,x_nbr)*hx' + 1/6*hxxx*reshape(dE_xf_xfxf_xf(:,:,jp2),x_nbr^3,x_nbr)*hx' + 1/6*hxxx*reshape(E_xf_xfxf_xf,x_nbr^3,x_nbr)*dhx(:,:,jp2)'... + 1/6*dhuuu(:,:,jp2)*reshape(QPu*E_u_u_u_u,u_nbr^3,u_nbr)*hu' + 1/6*huuu*reshape(dE_u_u_u_u_jp2,u_nbr^3,u_nbr)*hu' + 1/6*huuu*reshape(QPu*E_u_u_u_u,u_nbr^3,u_nbr)*dhu(:,:,jp2)'... diff --git a/matlab/+pruned_SS/quadruplication.m b/matlab/+pruned_SS/quadruplication.m index f6a31a29c..6350d9735 100644 --- a/matlab/+pruned_SS/quadruplication.m +++ b/matlab/+pruned_SS/quadruplication.m @@ -49,7 +49,7 @@ for l=1:p for k=l:p for j=k:p for i=j:p - idx = uperm([i j k l]); + idx = pruned_SS.uperm([i j k l]); for r = 1:size(idx,1) ii = idx(r,1); jj= idx(r,2); kk=idx(r,3); ll=idx(r,4); n = ii + (jj-1)*p + (kk-1)*p^2 + (ll-1)*p^3; diff --git a/matlab/uperm.m b/matlab/+pruned_SS/uperm.m similarity index 96% rename from matlab/uperm.m rename to matlab/+pruned_SS/uperm.m index 5b7e8df20..28a22d625 100644 --- a/matlab/uperm.m +++ b/matlab/+pruned_SS/uperm.m @@ -1,49 +1,49 @@ -function p = uperm(a) -% Return all unique permutations of possibly-repeating array elements -% ========================================================================= -% Copyright © 2014 Bruno Luong -% Copyright © 2020 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 . -% ========================================================================= -% Original author: Bruno Luong , April 20, 2014 -% https://groups.google.com/d/msg/comp.soft-sys.matlab/yQKVPTYrv6Q/gw1MzNd9sYkJ -% https://stackoverflow.com/a/42810388 - -[u, ~, J] = unique(a); -p = u(up(J, length(a))); - -function p = up(J, n) -ktab = histc(J,1:max(J)); -l = n; -p = zeros(1, n); -s = 1; -for i=1:length(ktab) - k = ktab(i); - c = nchoosek(1:l, k); - m = size(c,1); - [t, ~] = find(~p.'); - t = reshape(t, [], s); - c = t(c,:)'; - s = s*m; - r = repmat((1:s)',[1 k]); - q = accumarray([r(:) c(:)], i, [s n]); - p = repmat(p, [m 1]) + q; - l = l - k; -end -end - +function p = uperm(a) +% Return all unique permutations of possibly-repeating array elements +% ========================================================================= +% Copyright © 2014 Bruno Luong +% Copyright © 2020 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 . +% ========================================================================= +% Original author: Bruno Luong , April 20, 2014 +% https://groups.google.com/d/msg/comp.soft-sys.matlab/yQKVPTYrv6Q/gw1MzNd9sYkJ +% https://stackoverflow.com/a/42810388 + +[u, ~, J] = unique(a); +p = u(up(J, length(a))); + +function p = up(J, n) +ktab = histc(J,1:max(J)); +l = n; +p = zeros(1, n); +s = 1; +for i=1:length(ktab) + k = ktab(i); + c = nchoosek(1:l, k); + m = size(c,1); + [t, ~] = find(~p.'); + t = reshape(t, [], s); + c = t(c,:)'; + s = s*m; + r = repmat((1:s)',[1 k]); + q = accumarray([r(:) c(:)], i, [s n]); + p = repmat(p, [m 1]) + q; + l = l - k; +end +end + end % uperm \ No newline at end of file diff --git a/matlab/simul_static_model.m b/matlab/backward/simul_static_model.m similarity index 100% rename from matlab/simul_static_model.m rename to matlab/backward/simul_static_model.m diff --git a/matlab/cellofchar2mfile.m b/matlab/cellofchar2mfile.m deleted file mode 100644 index 0862058c1..000000000 --- a/matlab/cellofchar2mfile.m +++ /dev/null @@ -1,63 +0,0 @@ -function cellofchar2mfile(fname, c, cname) - -% Write a cell of char in a matlab script. -% -% INPUTS -% - fname [string] name of the file where c is to be saved. -% - c [cell] a two dimensional cell of char. -% -% OUTPUTS -% None. - -% Copyright © 2015-2017 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 . - -[pathstr,name,ext] = fileparts(fname); - -if isempty(ext) - fname = [pathstr, name, '.m'] -else - if ~isequal(ext, '.m') - error(['The first argument needs to be the name of a matlab script (with an .m extension)!']) - end -end - -if ~iscell(c) - error('The second input argument must be a cell!') -end - -if ndims(c)>2 - error(['The cell passed has a second argument cannot have more than two dimensions!']) -end - -variablename = inputname(2); - -if isempty(variablename) && nargin<3 - error(['You must pass the name of the cell (second input argument) as a string in the third input argument!']) -end - -if nargin>2 - if isvarname(cname) - variablename = cname; - else - error('The third input argument must be a valid variable name!') - end -end - -fid = fopen(fname,'w'); -fprintf(fid, '%s = %s;', variablename, writecellofchar(c)); -fclose(fid); \ No newline at end of file diff --git a/matlab/check_model.m b/matlab/check_model.m index 71dc6a2b5..7dc297109 100644 --- a/matlab/check_model.m +++ b/matlab/check_model.m @@ -2,7 +2,7 @@ function check_model(M_) % check_model(M_) % Performs various consistency checks on the model -% Copyright (C) 2005-2033 Dynare Team +% Copyright (C) 2005-2023 Dynare Team % % This file is part of Dynare. % diff --git a/matlab/clear_persistent_variables.m b/matlab/clear_persistent_variables.m index 989d4a361..32223192b 100644 --- a/matlab/clear_persistent_variables.m +++ b/matlab/clear_persistent_variables.m @@ -1,5 +1,5 @@ function clear_persistent_variables(folder, writelistofroutinestobecleared) - +% clear_persistent_variables(folder, writelistofroutinestobecleared) % Clear all the functions with persistent variables in directory folder (and subdirectories). % Copyright © 2015-2019 Dynare Team @@ -60,3 +60,51 @@ end list_of_functions_to_be_cleared; clear(list_of_functions{:}); + +function cellofchar2mfile(fname, c, cname) +% Write a cell of char in a matlab script. +% +% INPUTS +% - fname [string] name of the file where c is to be saved. +% - c [cell] a two dimensional cell of char. +% +% OUTPUTS +% None. + + + +[pathstr,name,ext] = fileparts(fname); + +if isempty(ext) + fname = [pathstr, name, '.m']; +else + if ~isequal(ext, '.m') + error(['The first argument needs to be the name of a matlab script (with an .m extension)!']) + end +end + +if ~iscell(c) + error('The second input argument must be a cell!') +end + +if ndims(c)>2 + error(['The cell passed has a second argument cannot have more than two dimensions!']) +end + +variablename = inputname(2); + +if isempty(variablename) && nargin<3 + error(['You must pass the name of the cell (second input argument) as a string in the third input argument!']) +end + +if nargin>2 + if isvarname(cname) + variablename = cname; + else + error('The third input argument must be a valid variable name!') + end +end + +fid = fopen(fname,'w'); +fprintf(fid, '%s = %s;', variablename, writecellofchar(c)); +fclose(fid); diff --git a/matlab/print_moments_implied_prior.m b/matlab/cli/print_moments_implied_prior.m similarity index 100% rename from matlab/print_moments_implied_prior.m rename to matlab/cli/print_moments_implied_prior.m diff --git a/matlab/compute_model_moments.m b/matlab/compute_model_moments.m deleted file mode 100644 index aa92954ff..000000000 --- a/matlab/compute_model_moments.m +++ /dev/null @@ -1,32 +0,0 @@ -function moments=compute_model_moments(dr,M_,options_) -% -% INPUTS -% dr: structure describing model solution -% M_: structure of Dynare options -% options_ -% -% OUTPUTS -% moments: a cell array containing requested moments -% -% SPECIAL REQUIREMENTS -% none - -% Copyright © 2008-2017 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 . - -[ivar,vartan,options_] = get_variables_list(options_,M_); -moments = th_autocovariances(dr,ivar,M_,options_,options_.nodecomposition); diff --git a/matlab/convergence_diagnostics/mcmc_ifac.m b/matlab/convergence_diagnostics/mcmc_ifac.m index 29d881dad..1c63bb7c2 100644 --- a/matlab/convergence_diagnostics/mcmc_ifac.m +++ b/matlab/convergence_diagnostics/mcmc_ifac.m @@ -69,3 +69,27 @@ for i=(Nc/2)+1: Nc+1 end Parzen=Parzen'; Ifac= 1+2*sum(Parzen(:).* AcorrXSIM); + + +function acf = dyn_autocorr(y, ar) +% function acf = dyn_autocorr(y, ar) +% autocorrelation function of y +% +% INPUTS +% y: time series +% ar: # of lags +% +% OUTPUTS +% acf: autocorrelation for lags 1 to ar +% +% SPECIAL REQUIREMENTS +% none + +y=y(:); +acf = NaN(ar+1,1); +acf(1)=1; +m = mean(y); +sd = std(y,1); +for i=1:ar + acf(i+1) = (y(i+1:end)-m)'*(y(1:end-i)-m)./((size(y,1))*sd^2); +end diff --git a/matlab/dyn_autocorr.m b/matlab/dyn_autocorr.m deleted file mode 100644 index dcf19ca31..000000000 --- a/matlab/dyn_autocorr.m +++ /dev/null @@ -1,40 +0,0 @@ -function acf = dyn_autocorr(y, ar) -% function acf = dyn_autocorr(y, ar) -% autocorrelation function of y -% -% INPUTS -% y: time series -% ar: # of lags -% -% OUTPUTS -% acf: autocorrelation for lags 1 to ar -% -% SPECIAL REQUIREMENTS -% none - -% Copyright © 2015-16 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 . - - -y=y(:); -acf = NaN(ar+1,1); -acf(1)=1; -m = mean(y); -sd = std(y,1); -for i=1:ar - acf(i+1) = (y(i+1:end)-m)'*(y(1:end-i)-m)./((size(y,1))*sd^2); -end diff --git a/matlab/dyn_diag_vech.m b/matlab/dyn_diag_vech.m deleted file mode 100644 index e07825cbb..000000000 --- a/matlab/dyn_diag_vech.m +++ /dev/null @@ -1,31 +0,0 @@ -function d = dyn_diag_vech(Vector) -% This function returns the diagonal elements of a symmetric matrix -% stored in vech form -% -% INPUTS -% Vector [double] a m*1 vector. -% -% OUTPUTS -% d [double] a n*1 vector, where n solves n*(n+1)/2=m. - -% Copyright © 2010-2017 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 . - -m = length(Vector); -n = (sqrt(1+8*m)-1)/2; -k = cumsum(1:n); -d = Vector(k); diff --git a/matlab/dynare_config.m b/matlab/dynare_config.m index 13f0d5c84..7702fd7eb 100644 --- a/matlab/dynare_config.m +++ b/matlab/dynare_config.m @@ -47,6 +47,7 @@ p = {'/../contrib/ms-sbvar/TZcode/MatlabFiles/' ; ... '/AIM/' ; ... '/backward/' ; ... '/cli/' ; ... + '/conditional_forecasts/'; ... '/convergence_diagnostics/' ; ... '/discretionary_policy/' ; ... '/distributions/' ; ... @@ -57,19 +58,24 @@ p = {'/../contrib/ms-sbvar/TZcode/MatlabFiles/' ; ... '/gsa/' ; ... '/kalman/' ; ... '/kalman/likelihood' ; ... + '/latex/' ; ... '/lmmcp/' ; ... '/modules/dseries/src/' ; ... '/reporting/' ; ... + '/matrix_solver/'; ... '/moments/'; ... '/ms-sbvar/' ; ... '/ms-sbvar/identification/' ; ... '/nonlinear-filters/' ; ... '/ols/' ; ... + '/optimal_policy/' ; ... '/optimization/' ; ... '/pac-tools/' ; ... '/parallel/' ; ... '/partial_information/' ; ... '/perfect-foresight-models/' ; ... + '/shock_decomposition/' ; ... + '/stochastic_solver/' ; ... '/utilities/dataset/' ; ... '/utilities/doc/' ; ... '/utilities/estimation/' ; ... diff --git a/matlab/dynare_gradient.m b/matlab/dynare_gradient.m deleted file mode 100644 index 641ad81d9..000000000 --- a/matlab/dynare_gradient.m +++ /dev/null @@ -1,66 +0,0 @@ -function [F,G] = dynare_gradient(fcn,x,epsilon,varargin) -% Computes the gradient of a function from R^m in R^n. -% -% INPUTS: -% fcn [string] name of the matlab's function. -% x [double] m*1 vector (where the gradient is evaluated). -% epsilon [double] scalar or m*1 vector of steps. -% -% OUTPUTS: -% F [double] n*1 vector, evaluation of the function at x. -% G [double] n*m matrix, evaluation of the gradient at x. -% -% OUTPUTS -% -% Copyright © 2010-2017 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 . - -% Evaluate the function at x. -F = feval(fcn, x, varargin{:}); - -% (G)Set dimensions. -m = length(x); -n = length(F); - -% Initialization of the gradient. -G = NaN(length(F),length(x)); - -if length(epsilon==1) - H = epsilon*eye(m); -else - H = diag(epsilon); -end - -% Compute the gradient. -for i=1:m - if size(x,1)>size(x,2) - h = H(i,:); - else - h = H(:,i); - end - [Fh,~,~,flag] = feval(fcn, x+transpose(h), varargin{:}); - if flag - G(:,i) = (Fh-F)/epsilon; - else - [Fh,~,~,flag] = feval(fcn, x-transpose(h), varargin{:}); - if flag - G(:,i) = (F-Fh)/epsilon; - else - error('-- Bad gradient --') - end - end -end \ No newline at end of file diff --git a/matlab/cartesian_product_of_sets.m b/matlab/ep/cartesian_product_of_sets.m similarity index 100% rename from matlab/cartesian_product_of_sets.m rename to matlab/ep/cartesian_product_of_sets.m diff --git a/matlab/cubature_with_gaussian_weight.m b/matlab/ep/cubature_with_gaussian_weight.m similarity index 100% rename from matlab/cubature_with_gaussian_weight.m rename to matlab/ep/cubature_with_gaussian_weight.m diff --git a/matlab/gauss_hermite_weights_and_nodes.m b/matlab/ep/gauss_hermite_weights_and_nodes.m similarity index 100% rename from matlab/gauss_hermite_weights_and_nodes.m rename to matlab/ep/gauss_hermite_weights_and_nodes.m diff --git a/matlab/gauss_legendre_weights_and_nodes.m b/matlab/ep/gauss_legendre_weights_and_nodes.m similarity index 100% rename from matlab/gauss_legendre_weights_and_nodes.m rename to matlab/ep/gauss_legendre_weights_and_nodes.m diff --git a/matlab/estimation/dsge_likelihood.m b/matlab/estimation/dsge_likelihood.m index 5ea603593..87689b61e 100644 --- a/matlab/estimation/dsge_likelihood.m +++ b/matlab/estimation/dsge_likelihood.m @@ -482,14 +482,14 @@ if analytic_derivation old_analytic_derivation_mode = options_.analytic_derivation_mode; options_.analytic_derivation_mode = kron_flag; if full_Hess - DERIVS = get_perturbation_params_derivs(M_, options_, estim_params_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, indparam, indexo, [], true); + DERIVS = identification.get_perturbation_params_derivs(M_, options_, estim_params_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, indparam, indexo, [], true); indD2T = reshape(1:M_.endo_nbr^2, M_.endo_nbr, M_.endo_nbr); indD2Om = dyn_unvech(1:M_.endo_nbr*(M_.endo_nbr+1)/2); D2T = DERIVS.d2KalmanA(indD2T(iv,iv),:); D2Om = DERIVS.d2Om(dyn_vech(indD2Om(iv,iv)),:); D2Yss = DERIVS.d2Yss(iv,:,:); else - DERIVS = get_perturbation_params_derivs(M_, options_, estim_params_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, indparam, indexo, [], false); + DERIVS = identification.get_perturbation_params_derivs(M_, options_, estim_params_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, indparam, indexo, [], false); end DT = zeros(M_.endo_nbr, M_.endo_nbr, size(DERIVS.dghx,3)); DT(:,M_.nstatic+(1:M_.nspred),:) = DERIVS.dghx; diff --git a/matlab/estimation/dsge_simulated_theoretical_covariance.m b/matlab/estimation/dsge_simulated_theoretical_covariance.m index 3fb978ada..93d9e1fef 100644 --- a/matlab/estimation/dsge_simulated_theoretical_covariance.m +++ b/matlab/estimation/dsge_simulated_theoretical_covariance.m @@ -132,7 +132,7 @@ for file = 1:NumberOfDrawsFiles if ~options_.pruning tmp = th_autocovariances(dr,ivar,M_,options_,nodecomposition); else - pruned_state_space = pruned_state_space_system(M_, options_, dr, obs_var, options_.ar, 1, 0); + pruned_state_space = pruned_SS.pruned_state_space_system(M_, options_, dr, obs_var, options_.ar, 1, 0); tmp{1} = pruned_state_space.Var_y; for i=1:nar tmp{i+1} = pruned_state_space.Corr_yi(:,:,i); diff --git a/matlab/evaluate_dynamic_model.m b/matlab/evaluate_dynamic_model.m deleted file mode 100644 index 6c57d6c1b..000000000 --- a/matlab/evaluate_dynamic_model.m +++ /dev/null @@ -1,29 +0,0 @@ -function residuals = evaluate_dynamic_model(dynamicmodel, endogenousvariables, exogenousvariables, params, steadystate, leadlagincidence, samplesize) - -% Copyright © 2016 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 . - -ny = length(steadystate); -periods = rows(exogenousvariables); - -residuals = zeros(ny,samplesize); -icols = find(leadlagincidence'); - -for t = 2:samplesize+1 - residuals(:,t-1) = dynamicmodel(endogenousvariables(icols), exogenousvariables, params, steadystate, t); - icols = icols + ny; -end \ No newline at end of file diff --git a/matlab/DsgeSmoother.m b/matlab/kalman/DsgeSmoother.m similarity index 100% rename from matlab/DsgeSmoother.m rename to matlab/kalman/DsgeSmoother.m diff --git a/matlab/evaluate_smoother.m b/matlab/kalman/evaluate_smoother.m similarity index 100% rename from matlab/evaluate_smoother.m rename to matlab/kalman/evaluate_smoother.m diff --git a/matlab/missing_DiffuseKalmanSmootherH1_Z.m b/matlab/kalman/missing_DiffuseKalmanSmootherH1_Z.m similarity index 100% rename from matlab/missing_DiffuseKalmanSmootherH1_Z.m rename to matlab/kalman/missing_DiffuseKalmanSmootherH1_Z.m diff --git a/matlab/missing_DiffuseKalmanSmootherH3_Z.m b/matlab/kalman/missing_DiffuseKalmanSmootherH3_Z.m similarity index 100% rename from matlab/missing_DiffuseKalmanSmootherH3_Z.m rename to matlab/kalman/missing_DiffuseKalmanSmootherH3_Z.m diff --git a/matlab/save_display_classical_smoother_results.m b/matlab/kalman/save_display_classical_smoother_results.m similarity index 100% rename from matlab/save_display_classical_smoother_results.m rename to matlab/kalman/save_display_classical_smoother_results.m diff --git a/matlab/store_smoother_results.m b/matlab/kalman/store_smoother_results.m similarity index 100% rename from matlab/store_smoother_results.m rename to matlab/kalman/store_smoother_results.m diff --git a/matlab/collect_latex_files.m b/matlab/latex/collect_latex_files.m similarity index 100% rename from matlab/collect_latex_files.m rename to matlab/latex/collect_latex_files.m diff --git a/matlab/dyn_latex_table.m b/matlab/latex/dyn_latex_table.m similarity index 100% rename from matlab/dyn_latex_table.m rename to matlab/latex/dyn_latex_table.m diff --git a/matlab/isbayes.m b/matlab/latex/isbayes.m similarity index 100% rename from matlab/isbayes.m rename to matlab/latex/isbayes.m diff --git a/matlab/write_latex_definitions.m b/matlab/latex/write_latex_definitions.m similarity index 100% rename from matlab/write_latex_definitions.m rename to matlab/latex/write_latex_definitions.m diff --git a/matlab/write_latex_parameter_table.m b/matlab/latex/write_latex_parameter_table.m similarity index 100% rename from matlab/write_latex_parameter_table.m rename to matlab/latex/write_latex_parameter_table.m diff --git a/matlab/write_latex_prior_table.m b/matlab/latex/write_latex_prior_table.m similarity index 100% rename from matlab/write_latex_prior_table.m rename to matlab/latex/write_latex_prior_table.m diff --git a/matlab/mcp_func.m b/matlab/lmmcp/mcp_func.m similarity index 100% rename from matlab/mcp_func.m rename to matlab/lmmcp/mcp_func.m diff --git a/matlab/fastgensylv.m b/matlab/matrix_solver/fastgensylv.m similarity index 100% rename from matlab/fastgensylv.m rename to matlab/matrix_solver/fastgensylv.m diff --git a/matlab/gensylv_fp.m b/matlab/matrix_solver/gensylv_fp.m similarity index 100% rename from matlab/gensylv_fp.m rename to matlab/matrix_solver/gensylv_fp.m diff --git a/matlab/lyapunov_solver.m b/matlab/matrix_solver/lyapunov_solver.m similarity index 100% rename from matlab/lyapunov_solver.m rename to matlab/matrix_solver/lyapunov_solver.m diff --git a/matlab/lyapunov_symm.m b/matlab/matrix_solver/lyapunov_symm.m similarity index 100% rename from matlab/lyapunov_symm.m rename to matlab/matrix_solver/lyapunov_symm.m diff --git a/matlab/quadratic_matrix_equation_solver.m b/matlab/matrix_solver/quadratic_matrix_equation_solver.m similarity index 100% rename from matlab/quadratic_matrix_equation_solver.m rename to matlab/matrix_solver/quadratic_matrix_equation_solver.m diff --git a/matlab/sylvester3.m b/matlab/matrix_solver/sylvester3.m similarity index 100% rename from matlab/sylvester3.m rename to matlab/matrix_solver/sylvester3.m diff --git a/matlab/sylvester3a.m b/matlab/matrix_solver/sylvester3a.m similarity index 100% rename from matlab/sylvester3a.m rename to matlab/matrix_solver/sylvester3a.m diff --git a/matlab/UnivariateSpectralDensity.m b/matlab/moments/UnivariateSpectralDensity.m similarity index 100% rename from matlab/UnivariateSpectralDensity.m rename to matlab/moments/UnivariateSpectralDensity.m diff --git a/matlab/add_filter_subtitle.m b/matlab/moments/add_filter_subtitle.m similarity index 100% rename from matlab/add_filter_subtitle.m rename to matlab/moments/add_filter_subtitle.m diff --git a/matlab/compute_moments_varendo.m b/matlab/moments/compute_moments_varendo.m similarity index 100% rename from matlab/compute_moments_varendo.m rename to matlab/moments/compute_moments_varendo.m diff --git a/matlab/conditional_variance_decomposition.m b/matlab/moments/conditional_variance_decomposition.m similarity index 100% rename from matlab/conditional_variance_decomposition.m rename to matlab/moments/conditional_variance_decomposition.m diff --git a/matlab/conditional_variance_decomposition_ME_mc_analysis.m b/matlab/moments/conditional_variance_decomposition_ME_mc_analysis.m similarity index 100% rename from matlab/conditional_variance_decomposition_ME_mc_analysis.m rename to matlab/moments/conditional_variance_decomposition_ME_mc_analysis.m diff --git a/matlab/conditional_variance_decomposition_mc_analysis.m b/matlab/moments/conditional_variance_decomposition_mc_analysis.m similarity index 100% rename from matlab/conditional_variance_decomposition_mc_analysis.m rename to matlab/moments/conditional_variance_decomposition_mc_analysis.m diff --git a/matlab/correlation_mc_analysis.m b/matlab/moments/correlation_mc_analysis.m similarity index 100% rename from matlab/correlation_mc_analysis.m rename to matlab/moments/correlation_mc_analysis.m diff --git a/matlab/covariance_mc_analysis.m b/matlab/moments/covariance_mc_analysis.m similarity index 100% rename from matlab/covariance_mc_analysis.m rename to matlab/moments/covariance_mc_analysis.m diff --git a/matlab/disp_moments.m b/matlab/moments/disp_moments.m similarity index 100% rename from matlab/disp_moments.m rename to matlab/moments/disp_moments.m diff --git a/matlab/moments/disp_th_moments_pruned_state_space.m b/matlab/moments/disp_th_moments_pruned_state_space.m index 8db2dad71..55d504458 100644 --- a/matlab/moments/disp_th_moments_pruned_state_space.m +++ b/matlab/moments/disp_th_moments_pruned_state_space.m @@ -52,7 +52,7 @@ for i=1:nvars obs_var(i,1) = find(strcmp(M_.endo_names(i_var(i),:), M_.endo_names(dr.order_var))); end -pruned_state_space = pruned_state_space_system(M_, options_, dr, obs_var, options_.ar, 1, 0); +pruned_state_space = pruned_SS.pruned_state_space_system(M_, options_, dr, obs_var, options_.ar, 1, 0); m = pruned_state_space.E_y; diff --git a/matlab/my_subplot.m b/matlab/my_subplot.m deleted file mode 100644 index ad03060a2..000000000 --- a/matlab/my_subplot.m +++ /dev/null @@ -1,56 +0,0 @@ -function my_subplot(i,imax,irow,icol,fig_title) - -% function my_subplot(i,imax,irow,icol,fig_title) -% spreads subplots on several figures according to a maximum number of -% subplots per figure -% -% INPUTS -% i: subplot number -% imax: total number of subplots -% irow: maximum number of rows in a figure -% icol: maximum number of columns in a figure -% fig_title: title to be repeated on each figure -% -% OUTPUT -% none -% -% SPECIAL REQUIREMENTS -% none - -% Copyright © 2003-2009 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 . - -nfig_max = irow*icol; -if imax < nfig_max - icol = ceil(sqrt(imax)); - irow=icol; - if (icol-1)*(icol-2) >= imax - irow = icol-2; - icol = icol-1; - elseif (icol)*(icol-2) >= imax - irow = icol-2; - elseif icol*(icol-1) >= imax - irow = icol-1; - end -end - -i1 = mod(i-1,nfig_max); -if i1 == 0 - figure('Name',fig_title); -end - -subplot(irow,icol,i1+1); \ No newline at end of file diff --git a/matlab/reduced_rank_cholesky.m b/matlab/nonlinear-filters/reduced_rank_cholesky.m similarity index 100% rename from matlab/reduced_rank_cholesky.m rename to matlab/nonlinear-filters/reduced_rank_cholesky.m diff --git a/matlab/dyn_ramsey_static.m b/matlab/optimal_policy/dyn_ramsey_static.m similarity index 100% rename from matlab/dyn_ramsey_static.m rename to matlab/optimal_policy/dyn_ramsey_static.m diff --git a/matlab/evaluate_planner_objective.m b/matlab/optimal_policy/evaluate_planner_objective.m similarity index 100% rename from matlab/evaluate_planner_objective.m rename to matlab/optimal_policy/evaluate_planner_objective.m diff --git a/matlab/get_optimal_policy_discount_factor.m b/matlab/optimal_policy/get_optimal_policy_discount_factor.m similarity index 100% rename from matlab/get_optimal_policy_discount_factor.m rename to matlab/optimal_policy/get_optimal_policy_discount_factor.m diff --git a/matlab/mult_elimination.m b/matlab/optimal_policy/mult_elimination.m similarity index 100% rename from matlab/mult_elimination.m rename to matlab/optimal_policy/mult_elimination.m diff --git a/matlab/csolve.m b/matlab/optimization/csolve.m similarity index 100% rename from matlab/csolve.m rename to matlab/optimization/csolve.m diff --git a/matlab/dynare_solve.m b/matlab/optimization/dynare_solve.m similarity index 100% rename from matlab/dynare_solve.m rename to matlab/optimization/dynare_solve.m diff --git a/matlab/lnsrch1.m b/matlab/optimization/lnsrch1.m similarity index 100% rename from matlab/lnsrch1.m rename to matlab/optimization/lnsrch1.m diff --git a/matlab/options2cell.m b/matlab/optimization/options2cell.m similarity index 100% rename from matlab/options2cell.m rename to matlab/optimization/options2cell.m diff --git a/matlab/solve1.m b/matlab/optimization/solve1.m similarity index 100% rename from matlab/solve1.m rename to matlab/optimization/solve1.m diff --git a/matlab/solve_one_boundary.m b/matlab/optimization/solve_one_boundary.m similarity index 100% rename from matlab/solve_one_boundary.m rename to matlab/optimization/solve_one_boundary.m diff --git a/matlab/step_length_correction.m b/matlab/optimization/step_length_correction.m similarity index 100% rename from matlab/step_length_correction.m rename to matlab/optimization/step_length_correction.m diff --git a/matlab/trust_region.m b/matlab/optimization/trust_region.m similarity index 100% rename from matlab/trust_region.m rename to matlab/optimization/trust_region.m diff --git a/matlab/basic_plan.m b/matlab/perfect-foresight-models/basic_plan.m similarity index 100% rename from matlab/basic_plan.m rename to matlab/perfect-foresight-models/basic_plan.m diff --git a/matlab/evaluate_max_dynamic_residual.m b/matlab/perfect-foresight-models/evaluate_max_dynamic_residual.m similarity index 100% rename from matlab/evaluate_max_dynamic_residual.m rename to matlab/perfect-foresight-models/evaluate_max_dynamic_residual.m diff --git a/matlab/flip_plan.m b/matlab/perfect-foresight-models/flip_plan.m similarity index 100% rename from matlab/flip_plan.m rename to matlab/perfect-foresight-models/flip_plan.m diff --git a/matlab/init_plan.m b/matlab/perfect-foresight-models/init_plan.m similarity index 100% rename from matlab/init_plan.m rename to matlab/perfect-foresight-models/init_plan.m diff --git a/matlab/save_results.m b/matlab/save_results.m deleted file mode 100644 index 7ee4801de..000000000 --- a/matlab/save_results.m +++ /dev/null @@ -1,38 +0,0 @@ -function save_results(x,s_name,names) - -% function save_results(x,s_name,names) -% save results in appropriate structure -% -% INPUT -% x: matrix to be saved column by column -% s_name: name of the structure where to save the results -% names: names of the individual series -% -% OUTPUT -% none -% -% SPECIAL REQUIREMENT -% none - -% Copyright © 2006-2009 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 . - -global oo_ - -for i=1:size(x,2) - eval([s_name deblank(names(i,:)) '= x(:,i);']); -end \ No newline at end of file diff --git a/matlab/WriteShockDecomp2Excel.m b/matlab/shock_decomposition/WriteShockDecomp2Excel.m similarity index 100% rename from matlab/WriteShockDecomp2Excel.m rename to matlab/shock_decomposition/WriteShockDecomp2Excel.m diff --git a/matlab/annualized_shock_decomposition.m b/matlab/shock_decomposition/annualized_shock_decomposition.m similarity index 100% rename from matlab/annualized_shock_decomposition.m rename to matlab/shock_decomposition/annualized_shock_decomposition.m diff --git a/matlab/epilogue_shock_decomposition.m b/matlab/shock_decomposition/epilogue_shock_decomposition.m similarity index 100% rename from matlab/epilogue_shock_decomposition.m rename to matlab/shock_decomposition/epilogue_shock_decomposition.m diff --git a/matlab/expand_group.m b/matlab/shock_decomposition/expand_group.m similarity index 100% rename from matlab/expand_group.m rename to matlab/shock_decomposition/expand_group.m diff --git a/matlab/graph_decomp.m b/matlab/shock_decomposition/graph_decomp.m similarity index 100% rename from matlab/graph_decomp.m rename to matlab/shock_decomposition/graph_decomp.m diff --git a/matlab/graph_decomp_detail.m b/matlab/shock_decomposition/graph_decomp_detail.m similarity index 100% rename from matlab/graph_decomp_detail.m rename to matlab/shock_decomposition/graph_decomp_detail.m diff --git a/matlab/initial_condition_decomposition.m b/matlab/shock_decomposition/initial_condition_decomposition.m similarity index 100% rename from matlab/initial_condition_decomposition.m rename to matlab/shock_decomposition/initial_condition_decomposition.m diff --git a/matlab/plot_shock_decomposition.m b/matlab/shock_decomposition/plot_shock_decomposition.m similarity index 100% rename from matlab/plot_shock_decomposition.m rename to matlab/shock_decomposition/plot_shock_decomposition.m diff --git a/matlab/realtime_shock_decomposition.m b/matlab/shock_decomposition/realtime_shock_decomposition.m similarity index 100% rename from matlab/realtime_shock_decomposition.m rename to matlab/shock_decomposition/realtime_shock_decomposition.m diff --git a/matlab/set_default_initial_condition_decomposition_options.m b/matlab/shock_decomposition/set_default_initial_condition_decomposition_options.m similarity index 100% rename from matlab/set_default_initial_condition_decomposition_options.m rename to matlab/shock_decomposition/set_default_initial_condition_decomposition_options.m diff --git a/matlab/set_default_plot_shock_decomposition_options.m b/matlab/shock_decomposition/set_default_plot_shock_decomposition_options.m similarity index 100% rename from matlab/set_default_plot_shock_decomposition_options.m rename to matlab/shock_decomposition/set_default_plot_shock_decomposition_options.m diff --git a/matlab/shock_decomposition.m b/matlab/shock_decomposition/shock_decomposition.m similarity index 100% rename from matlab/shock_decomposition.m rename to matlab/shock_decomposition/shock_decomposition.m diff --git a/matlab/squeeze_shock_decomposition.m b/matlab/shock_decomposition/squeeze_shock_decomposition.m similarity index 100% rename from matlab/squeeze_shock_decomposition.m rename to matlab/shock_decomposition/squeeze_shock_decomposition.m diff --git a/matlab/AIM_first_order_solver.m b/matlab/stochastic_solver/AIM_first_order_solver.m similarity index 100% rename from matlab/AIM_first_order_solver.m rename to matlab/stochastic_solver/AIM_first_order_solver.m diff --git a/matlab/compute_decision_rules.m b/matlab/stochastic_solver/compute_decision_rules.m similarity index 100% rename from matlab/compute_decision_rules.m rename to matlab/stochastic_solver/compute_decision_rules.m diff --git a/matlab/convertAimCodeToInfo.m b/matlab/stochastic_solver/convertAimCodeToInfo.m similarity index 100% rename from matlab/convertAimCodeToInfo.m rename to matlab/stochastic_solver/convertAimCodeToInfo.m diff --git a/matlab/disp_dr.m b/matlab/stochastic_solver/disp_dr.m similarity index 100% rename from matlab/disp_dr.m rename to matlab/stochastic_solver/disp_dr.m diff --git a/matlab/dyn_first_order_solver.m b/matlab/stochastic_solver/dyn_first_order_solver.m similarity index 100% rename from matlab/dyn_first_order_solver.m rename to matlab/stochastic_solver/dyn_first_order_solver.m diff --git a/matlab/dyn_second_order_solver.m b/matlab/stochastic_solver/dyn_second_order_solver.m similarity index 100% rename from matlab/dyn_second_order_solver.m rename to matlab/stochastic_solver/dyn_second_order_solver.m diff --git a/matlab/dynare_resolve.m b/matlab/stochastic_solver/dynare_resolve.m similarity index 100% rename from matlab/dynare_resolve.m rename to matlab/stochastic_solver/dynare_resolve.m diff --git a/matlab/getIrfShocksIndx.m b/matlab/stochastic_solver/getIrfShocksIndx.m similarity index 100% rename from matlab/getIrfShocksIndx.m rename to matlab/stochastic_solver/getIrfShocksIndx.m diff --git a/matlab/irf.m b/matlab/stochastic_solver/irf.m similarity index 100% rename from matlab/irf.m rename to matlab/stochastic_solver/irf.m diff --git a/matlab/k_order_pert.m b/matlab/stochastic_solver/k_order_pert.m similarity index 100% rename from matlab/k_order_pert.m rename to matlab/stochastic_solver/k_order_pert.m diff --git a/matlab/kalman_transition_matrix.m b/matlab/stochastic_solver/kalman_transition_matrix.m similarity index 100% rename from matlab/kalman_transition_matrix.m rename to matlab/stochastic_solver/kalman_transition_matrix.m diff --git a/matlab/resol.m b/matlab/stochastic_solver/resol.m similarity index 100% rename from matlab/resol.m rename to matlab/stochastic_solver/resol.m diff --git a/matlab/select_qz_criterium_value.m b/matlab/stochastic_solver/select_qz_criterium_value.m similarity index 100% rename from matlab/select_qz_criterium_value.m rename to matlab/stochastic_solver/select_qz_criterium_value.m diff --git a/matlab/set_state_space.m b/matlab/stochastic_solver/set_state_space.m similarity index 100% rename from matlab/set_state_space.m rename to matlab/stochastic_solver/set_state_space.m diff --git a/matlab/simult.m b/matlab/stochastic_solver/simult.m similarity index 100% rename from matlab/simult.m rename to matlab/stochastic_solver/simult.m diff --git a/matlab/simult_.m b/matlab/stochastic_solver/simult_.m similarity index 100% rename from matlab/simult_.m rename to matlab/stochastic_solver/simult_.m diff --git a/matlab/simultxdet.m b/matlab/stochastic_solver/simultxdet.m similarity index 100% rename from matlab/simultxdet.m rename to matlab/stochastic_solver/simultxdet.m diff --git a/matlab/stoch_simul.m b/matlab/stochastic_solver/stoch_simul.m similarity index 100% rename from matlab/stoch_simul.m rename to matlab/stochastic_solver/stoch_simul.m diff --git a/matlab/stochastic_solvers.m b/matlab/stochastic_solver/stochastic_solvers.m similarity index 100% rename from matlab/stochastic_solvers.m rename to matlab/stochastic_solver/stochastic_solvers.m diff --git a/tests/analytic_derivatives/BrockMirman_PertParamsDerivs.mod b/tests/analytic_derivatives/BrockMirman_PertParamsDerivs.mod index 0c5ac994d..9c181990a 100644 --- a/tests/analytic_derivatives/BrockMirman_PertParamsDerivs.mod +++ b/tests/analytic_derivatives/BrockMirman_PertParamsDerivs.mod @@ -440,7 +440,7 @@ for jj = 1:2 oo_.dr.Correlation_matrix = M_.Correlation_matrix; ex0 = oo_.exo_steady_state'; [~, oo_.dr.g1, oo_.dr.g2, oo_.dr.g3] = feval([M_.fname,'.dynamic'], oo_.dr.ys(I), oo_.exo_steady_state', M_.params, oo_.dr.ys, 1); - oo_.dr.g3 = unfold_g3(oo_.dr.g3, length(oo_.dr.ys(I))+length(oo_.exo_steady_state')); %add symmetric elements to g3 + oo_.dr.g3 = identification.unfold_g3(oo_.dr.g3, length(oo_.dr.ys(I))+length(oo_.exo_steady_state')); %add symmetric elements to g3 fprintf('***** %s: SOME COMMON OBJECTS *****\n', strparamset) for id_var = 1:size(lst_vars,2) @@ -460,7 +460,7 @@ for jj = 1:2 for id_kronflag = 1:length(KRONFLAG) fprintf('***** %s: d2flag=%d and kronflag=%d *****\n',strparamset, d2flag,KRONFLAG(id_kronflag)) options_.analytic_derivation_mode = KRONFLAG(id_kronflag); - DERIVS = get_perturbation_params_derivs(M_, options_, estim_params_, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state, indpmodel, indpstderr, indpcorr, d2flag); + DERIVS = identification.get_perturbation_params_derivs(M_, options_, estim_params_, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state, indpmodel, indpstderr, indpcorr, d2flag); for id_var = 1:size(lst_dvars,2) dx = norm( vec(nSYM.(sprintf('%s',lst_dvars{id_var}))) - vec(DERIVS.(sprintf('%s',lst_dvars{id_var}))), Inf); fprintf('Max absolute deviation for %s: %e\n', lst_dvars{id_var}, dx); diff --git a/tests/analytic_derivatives/burnside_3_order_PertParamsDerivs.mod b/tests/analytic_derivatives/burnside_3_order_PertParamsDerivs.mod index ba7f17762..b1f3351ed 100644 --- a/tests/analytic_derivatives/burnside_3_order_PertParamsDerivs.mod +++ b/tests/analytic_derivatives/burnside_3_order_PertParamsDerivs.mod @@ -170,7 +170,7 @@ KRONFLAGS = [-1 -2 0 1]; for k = 1:length(KRONFLAGS) fprintf('KRONFLAG=%d\n',KRONFLAGS(k)); options_.analytic_derivation_mode = KRONFLAGS(k); - DERIVS = get_perturbation_params_derivs(M_, options_, estim_params_, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state, indpmodel, indpstderr, indpcorr, d2flag); + DERIVS = identification.get_perturbation_params_derivs(M_, options_, estim_params_, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state, indpmodel, indpstderr, indpcorr, d2flag); oo_.dr.dg_0 = permute(1/2*DERIVS.dghs2,[1 3 2]); oo_.dr.dg_1 = cat(2,DERIVS.dghx,DERIVS.dghu) + 3/6*cat(2,DERIVS.dghxss,DERIVS.dghuss); oo_.dr.dg_2 = 1/2*cat(2,DERIVS.dghxx,DERIVS.dghxu,DERIVS.dghuu); diff --git a/tests/minimal_state_space_system/as2007_minimal.mod b/tests/minimal_state_space_system/as2007_minimal.mod index 81bb6c049..09b924615 100644 --- a/tests/minimal_state_space_system/as2007_minimal.mod +++ b/tests/minimal_state_space_system/as2007_minimal.mod @@ -81,7 +81,7 @@ SS.B = oo_.dr.ghu(indx,:); SS.C = oo_.dr.ghx(indy,:); SS.D = oo_.dr.ghu(indy,:); -[CheckCO,minnx,minSS] = get_minimal_state_representation(SS,0); +[CheckCO,minnx,minSS] = identification.get_minimal_state_representation(SS,0); Sigmax_full = lyapunov_symm(SS.A, SS.B*M_.Sigma_e*SS.B', options_.lyapunov_fixed_point_tol, options_.qz_criterium, options_.lyapunov_complex_threshold, 1, options_.debug); Sigmay_full = SS.C*Sigmax_full*SS.C' + SS.D*M_.Sigma_e*SS.D'; diff --git a/tests/minimal_state_space_system/sw_minimal.mod b/tests/minimal_state_space_system/sw_minimal.mod index 585e851f3..b7dc1150b 100644 --- a/tests/minimal_state_space_system/sw_minimal.mod +++ b/tests/minimal_state_space_system/sw_minimal.mod @@ -414,7 +414,7 @@ SS.B = oo_.dr.ghu(indx,:); SS.C = oo_.dr.ghx(indy,:); SS.D = oo_.dr.ghu(indy,:); -[CheckCO,minnx,minSS] = get_minimal_state_representation(SS,0); +[CheckCO,minnx,minSS] = identification.get_minimal_state_representation(SS,0); Sigmax_full = lyapunov_symm(SS.A, SS.B*M_.Sigma_e*SS.B', options_.lyapunov_fixed_point_tol, options_.qz_criterium, options_.lyapunov_complex_threshold, 1, options_.debug); Sigmay_full = SS.C*Sigmax_full*SS.C' + SS.D*M_.Sigma_e*SS.D'; diff --git a/tests/stochastic_simulations/pruning/AS_pruned_state_space_red_shock.mod b/tests/stochastic_simulations/pruning/AS_pruned_state_space_red_shock.mod index 08c96ecdb..e2dc1af98 100644 --- a/tests/stochastic_simulations/pruning/AS_pruned_state_space_red_shock.mod +++ b/tests/stochastic_simulations/pruning/AS_pruned_state_space_red_shock.mod @@ -97,7 +97,7 @@ steady; check; model_diagnostics; @#for orderApp in [1, 2, 3] stoch_simul(order=@{orderApp},pruning,irf=0,periods=0); - pruned_state_space.order_@{orderApp} = pruned_state_space_system(M_, options_, oo_.dr, [], options_.ar, 1, 0); + pruned_state_space.order_@{orderApp} = pruned_SS.pruned_state_space_system(M_, options_, oo_.dr, [], options_.ar, 1, 0); @#if Andreasen_et_al_toolbox addpath('Dynare44Pruning_v2/simAndMoments3order'); %provide path to toolbox optPruning.orderApp = @{orderApp}; diff --git a/tests/stochastic_simulations/pruning/AnSchorfheide_pruned_state_space.mod b/tests/stochastic_simulations/pruning/AnSchorfheide_pruned_state_space.mod index 623b966cf..39d644ac7 100644 --- a/tests/stochastic_simulations/pruning/AnSchorfheide_pruned_state_space.mod +++ b/tests/stochastic_simulations/pruning/AnSchorfheide_pruned_state_space.mod @@ -96,7 +96,7 @@ steady; check; model_diagnostics; @#for orderApp in [1, 2, 3] stoch_simul(order=@{orderApp},pruning,irf=0,periods=0); - pruned_state_space.order_@{orderApp} = pruned_state_space_system(M_, options_, oo_.dr, [], options_.ar, 1, 0); + pruned_state_space.order_@{orderApp} = pruned_SS.pruned_state_space_system(M_, options_, oo_.dr, [], options_.ar, 1, 0); @#if Andreasen_et_al_toolbox addpath('Dynare44Pruning_v2/simAndMoments3order'); %provide path to toolbox optPruning.orderApp = @{orderApp}; diff --git a/tests/stochastic_simulations/pruning/fs2000_pruning.mod b/tests/stochastic_simulations/pruning/fs2000_pruning.mod index f325bf71f..732e4be18 100644 --- a/tests/stochastic_simulations/pruning/fs2000_pruning.mod +++ b/tests/stochastic_simulations/pruning/fs2000_pruning.mod @@ -97,7 +97,7 @@ if (max(abs(Y2_local(:) - Y2_simult(:)))>1e-12) end % pruned_state_space_system.m: implements Andreasen et al. -pss = pruned_state_space_system(M_, options_, oo_.dr, [], 0, false, false); +pss = pruned_SS.pruned_state_space_system(M_, options_, oo_.dr, [], 0, false, false); Y2_an = zeros(M_.endo_nbr,T+1); Y2_an(:,1) = oo_.dr.ys; % z = [xf;xs;kron(xf,xf)] @@ -127,7 +127,7 @@ stoch_simul(order=3, nograph, irf=0); % simult_.m Y3_simult = simult_(M_,options_,oo_.dr.ys,oo_.dr,ex,options_.order); % pruned_state_space_system.m -pss = pruned_state_space_system(M_, options_, oo_.dr, [], 0, false, false); +pss = pruned_SS.pruned_state_space_system(M_, options_, oo_.dr, [], 0, false, false); Y3_an = zeros(M_.endo_nbr,T+1); Y3_an(:,1) = oo_.dr.ys; % z = [xf; xs; kron(xf,xf); xrd; kron(xf,xs); kron(kron(xf,xf),xf)] From f05a2de89ee10bf14fe57308fce11dea65b7e9d9 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Mon, 11 Dec 2023 11:48:34 +0100 Subject: [PATCH 3/3] get_perturbation_params_derivs.m: replace try-catch by proper check of file existence Let's other errors though with explicit message --- .../get_perturbation_params_derivs.m | 66 +++++++------------ 1 file changed, 24 insertions(+), 42 deletions(-) diff --git a/matlab/+identification/get_perturbation_params_derivs.m b/matlab/+identification/get_perturbation_params_derivs.m index 9f6543998..579ebce38 100644 --- a/matlab/+identification/get_perturbation_params_derivs.m +++ b/matlab/+identification/get_perturbation_params_derivs.m @@ -447,13 +447,15 @@ if analytic_derivation_mode == -2 clear dYss_g elseif (analytic_derivation_mode == 0 || analytic_derivation_mode == 1) - %% Analytical computation of Jacobian and Hessian (wrt selected model parameters) of steady state, i.e. dYss and d2Yss - [~, g1_static] = feval([fname,'.static'], ys, exo_steady_state', params); %g1_static is [endo_nbr by endo_nbr] first-derivative (wrt all endogenous variables) of static model equations f, i.e. df/dys, in declaration order - try - rp_static = feval([fname,'.static_params_derivs'], ys, exo_steady_state', params); %rp_static is [endo_nbr by param_nbr] first-derivative (wrt all model parameters) of static model equations f, i.e. df/dparams, in declaration order - catch + if ~exist(['+' fname filesep 'static_params_derivs.m'],'file') error('For analytical parameter derivatives ''static_params_derivs.m'' file is needed, this can be created by putting identification(order=%d) into your mod file.',order) end + if ~exist(['+' fname filesep 'dynamic_params_derivs.m'],'file') + error('For analytical parameter derivatives ''dynamic_params_derivs.m'' file is needed, this can be created by putting identification(order=%d) into your mod file.',order) + end + %% Analytical computation of Jacobian and Hessian (wrt selected model parameters) of steady state, i.e. dYss and d2Yss + [~, g1_static] = feval([fname,'.static'], ys, exo_steady_state', params); %g1_static is [endo_nbr by endo_nbr] first-derivative (wrt all endogenous variables) of static model equations f, i.e. df/dys, in declaration order + rp_static = feval([fname,'.static_params_derivs'], ys, exo_steady_state', params); %rp_static is [endo_nbr by param_nbr] first-derivative (wrt all model parameters) of static model equations f, i.e. df/dparams, in declaration order dys = -g1_static\rp_static; %use implicit function theorem (equation 5 of Ratto and Iskrev (2012) to compute [endo_nbr by param_nbr] first-derivative (wrt all model parameters) of steady state for all endogenous variables analytically, note that dys is in declaration order d2ys = zeros(endo_nbr, param_nbr, param_nbr); %initialize in tensor notation, note that d2ys is only needed for d2flag, i.e. for g1pp if d2flag @@ -474,14 +476,10 @@ elseif (analytic_derivation_mode == 0 || analytic_derivation_mode == 1) %g1 is [endo_nbr by yy0ex0_nbr first derivative (wrt all dynamic variables) of dynamic model equations, i.e. df/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order %g2 is [endo_nbr by yy0ex0_nbr^2] second derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order %g3 is [endo_nbr by yy0ex0_nbr^3] third-derivative (wrt all dynamic variables) of dynamic model equations, i.e. (d(df/dyy0ex0)/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order - try - [~, g1p_static, rpp_static] = feval([fname,'.static_params_derivs'], ys, exo_steady_state', params); - %g1p_static is [endo_nbr by endo_nbr by param_nbr] first derivative (wrt all model parameters) of first-derivative (wrt all endogenous variables) of static model equations f, i.e. (df/dys)/dparams, in declaration order - %rpp_static is [#second_order_residual_terms by 4] and contains nonzero values and corresponding indices of second derivatives (wrt all model parameters) of static model equations f, i.e. d(df/dparams)/dparams, in declaration order, where - % column 1 contains equation number; column 2 contains first parameter; column 3 contains second parameter; column 4 contains value of derivative - catch - error('For analytical parameter derivatives ''static_params_derivs.m'' file is needed, this can be created by putting identification(order=%d) into your mod file.',order) - end + [~, g1p_static, rpp_static] = feval([fname,'.static_params_derivs'], ys, exo_steady_state', params); + %g1p_static is [endo_nbr by endo_nbr by param_nbr] first derivative (wrt all model parameters) of first-derivative (wrt all endogenous variables) of static model equations f, i.e. (df/dys)/dparams, in declaration order + %rpp_static is [#second_order_residual_terms by 4] and contains nonzero values and corresponding indices of second derivatives (wrt all model parameters) of static model equations f, i.e. d(df/dparams)/dparams, in declaration order, where + % column 1 contains equation number; column 2 contains first parameter; column 3 contains second parameter; column 4 contains value of derivative rpp_static = get_all_resid_2nd_derivs(rpp_static, endo_nbr, param_nbr); %make full matrix out of nonzero values and corresponding indices %rpp_static is [endo_nbr by param_nbr by param_nbr] second derivatives (wrt all model parameters) of static model equations, i.e. d(df/dparams)/dparams, in declaration order if isempty(find(g2_static)) @@ -525,50 +523,34 @@ elseif (analytic_derivation_mode == 0 || analytic_derivation_mode == 1) end if d2flag - try - if order < 3 - [~, g1p, ~, g1pp, g2p] = feval([fname,'.dynamic_params_derivs'], ys(I), exo_steady_state', params, ys, 1, dys, d2ys); - else - [~, g1p, ~, g1pp, g2p, g3p] = feval([fname,'.dynamic_params_derivs'], ys(I), exo_steady_state', params, ys, 1, dys, d2ys); - end - catch - error('For analytical parameter derivatives ''dynamic_params_derivs.m'' file is needed, this can be created by putting identification(order=%d) into your mod file.',order) + if order < 3 + [~, g1p, ~, g1pp, g2p] = feval([fname,'.dynamic_params_derivs'], ys(I), exo_steady_state', params, ys, 1, dys, d2ys); + else + [~, g1p, ~, g1pp, g2p, g3p] = feval([fname,'.dynamic_params_derivs'], ys(I), exo_steady_state', params, ys, 1, dys, d2ys); end %g1pp are nonzero values and corresponding indices of second-derivatives (wrt all model parameters) of first-derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(d(df/dyy0ex0)/dparam)/dparam, rows are in declaration order, first column in declaration order d2Yss = d2ys(order_var,indpmodel,indpmodel); %[endo_nbr by mod_param_nbr by mod_param_nbr], put into DR order and focus only on selected model parameters else if order == 1 - try - [~, g1p] = feval([fname,'.dynamic_params_derivs'], ys(I), exo_steady_state', params, ys, 1, dys, d2ys); - %g1p is [endo_nbr by yy0ex0_nbr by param_nbr] first-derivative (wrt all model parameters) of first-derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dparam, rows are in declaration order, column in lead_lag_incidence order - catch - error('For analytical parameter derivatives ''dynamic_params_derivs.m'' file is needed, this can be created by putting identification(order=%d) into your mod file.',order) - end + [~, g1p] = feval([fname,'.dynamic_params_derivs'], ys(I), exo_steady_state', params, ys, 1, dys, d2ys); + %g1p is [endo_nbr by yy0ex0_nbr by param_nbr] first-derivative (wrt all model parameters) of first-derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dparam, rows are in declaration order, column in lead_lag_incidence order [~, g1, g2 ] = feval([fname,'.dynamic'], ys(I), exo_steady_state', params, ys, 1); %g1 is [endo_nbr by yy0ex0_nbr first derivative (wrt all dynamic variables) of dynamic model equations, i.e. df/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order %g2 is [endo_nbr by yy0ex0_nbr^2] second derivatives (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order elseif order == 2 - try - [~, g1p, ~, ~, g2p] = feval([fname,'.dynamic_params_derivs'], ys(I), exo_steady_state', params, ys, 1, dys, d2ys); - %g1p is [endo_nbr by yy0ex0_nbr by param_nbr] first-derivative (wrt all model parameters) of first-derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dparam, rows are in declaration order, column in lead_lag_incidence order - %g2p are nonzero values and corresponding indices of first-derivative (wrt all model parameters) of second-derivatives (wrt all dynamic variables) of dynamic model equations, i.e. d(d(df/dyy0ex0)/dyy0ex0)/dparam, rows are in declaration order, first and second column in declaration order - catch - error('For analytical parameter derivatives ''dynamic_params_derivs.m'' file is needed, this can be created by putting identification(order=%d) into your mod file.',order) - end + [~, g1p, ~, ~, g2p] = feval([fname,'.dynamic_params_derivs'], ys(I), exo_steady_state', params, ys, 1, dys, d2ys); + %g1p is [endo_nbr by yy0ex0_nbr by param_nbr] first-derivative (wrt all model parameters) of first-derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dparam, rows are in declaration order, column in lead_lag_incidence order + %g2p are nonzero values and corresponding indices of first-derivative (wrt all model parameters) of second-derivatives (wrt all dynamic variables) of dynamic model equations, i.e. d(d(df/dyy0ex0)/dyy0ex0)/dparam, rows are in declaration order, first and second column in declaration order [~, g1, g2, g3] = feval([fname,'.dynamic'], ys(I), exo_steady_state', params, ys, 1); %note that g3 does not contain symmetric elements g3 = identification.unfold_g3(g3, yy0ex0_nbr); %add symmetric elements to g3 %g1 is [endo_nbr by yy0ex0_nbr first derivative (wrt all dynamic variables) of dynamic model equations, i.e. df/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order %g2 is [endo_nbr by yy0ex0_nbr^2] second derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order %g3 is [endo_nbr by yy0ex0_nbr^3] third-derivative (wrt all dynamic variables) of dynamic model equations, i.e. (d(df/dyy0ex0)/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order elseif order == 3 - try - [~, g1p, ~, ~, g2p, g3p] = feval([fname,'.dynamic_params_derivs'], ys(I), exo_steady_state', params, ys, 1, dys, d2ys); - %g1p is [endo_nbr by yy0ex0_nbr by param_nbr] first-derivative (wrt all model parameters) of first-derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dparam, rows are in declaration order, column in lead_lag_incidence order - %g2p are nonzero values and corresponding indices of first-derivative (wrt all model parameters) of second-derivatives (wrt all dynamic variables) of dynamic model equations, i.e. d(d(df/dyy0ex0)/dyy0ex0)/dparam, rows are in declaration order, first and second column in declaration order - %g3p are nonzero values and corresponding indices of first-derivative (wrt all model parameters) of third-derivatives (wrt all dynamic variables) of dynamic model equations, i.e. d(d(d(df/dyy0ex0)/dyy0ex0)/dyy0ex0)/dparam, rows are in declaration order, first, second and third column in declaration order - catch - error('For analytical parameter derivatives ''dynamic_params_derivs.m'' file is needed, this can be created by putting identification(order=%d) into your mod file.',order) - end + [~, g1p, ~, ~, g2p, g3p] = feval([fname,'.dynamic_params_derivs'], ys(I), exo_steady_state', params, ys, 1, dys, d2ys); + %g1p is [endo_nbr by yy0ex0_nbr by param_nbr] first-derivative (wrt all model parameters) of first-derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dparam, rows are in declaration order, column in lead_lag_incidence order + %g2p are nonzero values and corresponding indices of first-derivative (wrt all model parameters) of second-derivatives (wrt all dynamic variables) of dynamic model equations, i.e. d(d(df/dyy0ex0)/dyy0ex0)/dparam, rows are in declaration order, first and second column in declaration order + %g3p are nonzero values and corresponding indices of first-derivative (wrt all model parameters) of third-derivatives (wrt all dynamic variables) of dynamic model equations, i.e. d(d(d(df/dyy0ex0)/dyy0ex0)/dyy0ex0)/dparam, rows are in declaration order, first, second and third column in declaration order T = NaN(sum(dynamic_tmp_nbr(1:5))); T = feval([fname, '.dynamic_g4_tt'], T, ys(I), exo_steady_state', params, ys, 1); g1 = feval([fname, '.dynamic_g1'], T, ys(I), exo_steady_state', params, ys, 1, false); %g1 is [endo_nbr by yy0ex0_nbr first derivative (wrt all dynamic variables) of dynamic model equations, i.e. df/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order