stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge branch 'posterior_moments' of git.dynare.org:JohannesPfeifer/dynare 

Ref. !2242
covariance-quadratic-approximation
Sébastien Villemot 2023-12-20 14:03:18 +01:00
parent fb18777e58 774cb18aac
commit 38d1c0538a
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
14 changed files with 55 additions and 55 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
matlab/+gsa/map_calibration.m
View File

@ -212,7 +212,7 @@ if ~isempty(indx_irf)
indx_irf_matrix(:,plot_indx(ij)) = indx_irf_matrix(:,plot_indx(ij)) + indx_irf(:,ij);
for ik=1:size(mat_irf{ij},2)
[Mean,Median,Var,HPD,Distrib] = ...
posterior_moments(mat_irf{ij}(:,ik),0,options_.mh_conf_sig);
posterior_moments(mat_irf{ij}(:,ik),options_.mh_conf_sig);
irf_mean{plot_indx(ij)} = [irf_mean{plot_indx(ij)}; Mean];
irf_median{plot_indx(ij)} = [irf_median{plot_indx(ij)}; Median];
irf_var{plot_indx(ij)} = [irf_var{plot_indx(ij)}; Var];
@ -417,7 +417,7 @@ if ~isempty(indx_moment)
indx_moment_matrix(:,plot_indx(ij)) = indx_moment_matrix(:,plot_indx(ij)) + indx_moment(:,ij);
for ik=1:size(mat_moment{ij},2)
[Mean,Median,Var,HPD,Distrib] = ...
posterior_moments(mat_moment{ij}(:,ik),0,options_.mh_conf_sig);
posterior_moments(mat_moment{ij}(:,ik),options_.mh_conf_sig);
moment_mean{plot_indx(ij)} = [moment_mean{plot_indx(ij)}; Mean];
moment_median{plot_indx(ij)} = [moment_median{plot_indx(ij)}; Median];
moment_var{plot_indx(ij)} = [moment_var{plot_indx(ij)}; Var];

2
matlab/+gsa/reduced_form_mapping.m
View File

@ -720,7 +720,7 @@ function indmcf = redform_mcf(y0, x0, options_mcf, options_, fname, parnames, es
hh_fig=dyn_figure(options_.nodisplay,'name',options_mcf.amcf_title);
[~, ~, ~, ~, post_deciles] = posterior_moments(y0,1,0.9);
[~, ~, ~, ~, post_deciles] = posterior_moments(y0,0.9);
post_deciles = [-inf; post_deciles; inf];
for jt=1:10

2
matlab/+identification/plot.m
View File

@ -493,7 +493,7 @@ else
post_median=NaN(1,nparam);
hpd_interval=NaN(nparam,2);
for i=1:nparam
[~, post_median(:,i), ~, hpd_interval(i,:)] = posterior_moments(VVV(:,i),0,0.9);
[~, post_median(:,i), ~, hpd_interval(i,:)] = posterior_moments(VVV(:,i),0.9);
end
bar(post_median)
hold on, plot(hpd_interval,'--*r'),

4
matlab/+occbin/forecast.m
View File

@ -88,7 +88,7 @@ if opts.replic
for i=1:M_.endo_nbr
for j=1:forecast
[post_mean(j,1), post_median(j,1), post_var(j,1), hpd_interval(j,:), post_deciles(j,:)] = posterior_moments(squeeze(zlin0(j,i,:)),1,options_.forecasts.conf_sig);
[post_mean(j,1), post_median(j,1), post_var(j,1), hpd_interval(j,:), post_deciles(j,:)] = posterior_moments(squeeze(zlin0(j,i,:)),options_.forecasts.conf_sig);
end
oo_.occbin.linear_forecast.Mean.(M_.endo_names{i})=post_mean;
oo_.occbin.linear_forecast.Median.(M_.endo_names{i})=post_median;
@ -99,7 +99,7 @@ if opts.replic
oo_.occbin.linear_forecast.Min.(M_.endo_names{i})=zlinmin(:,i);
oo_.occbin.linear_forecast.Max.(M_.endo_names{i})=zlinmax(:,i);
for j=1:forecast
[post_mean(j,1), post_median(j,1), post_var(j,1), hpd_interval(j,:), post_deciles(j,:)] = posterior_moments(squeeze(zpiece0(j,i,:)),1,options_.forecasts.conf_sig);
[post_mean(j,1), post_median(j,1), post_var(j,1), hpd_interval(j,:), post_deciles(j,:)] = posterior_moments(squeeze(zpiece0(j,i,:)),options_.forecasts.conf_sig);
end
oo_.occbin.forecast.Mean.(M_.endo_names{i})=post_mean;
oo_.occbin.forecast.Median.(M_.endo_names{i})=post_median;

20
matlab/estimation/GetPosteriorParametersStatistics.m
View File

@ -105,7 +105,7 @@ if np
for i=1:np
if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh) || ishssmc(options_)
draws = getalldraws(ip);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, 1, options_.mh_conf_sig);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, options_.mh_conf_sig);
name = bayestopt_.name{ip};
oo_ = Filloo(oo_, name, type, post_mean, hpd_interval, post_median, post_var, post_deciles, density);
else
@ -114,7 +114,7 @@ if np
[post_mean, hpd_interval, post_var] = Extractoo(oo_, name, type);
catch
draws = getalldraws(ip);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, 1, options_.mh_conf_sig);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, options_.mh_conf_sig);
name = bayestopt_.name{ip};
oo_ = Filloo(oo_, name, type, post_mean, hpd_interval, post_median, post_var, post_deciles, density);
end
@ -148,7 +148,7 @@ if nvx
for i=1:nvx
if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh)
draws = getalldraws(ip);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, 1, options_.mh_conf_sig);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, options_.mh_conf_sig);
k = estim_params_.var_exo(i,1);
name = M_.exo_names{k};
oo_ = Filloo(oo_, name, type, post_mean, hpd_interval, post_median, post_var, post_deciles, density);
@ -160,7 +160,7 @@ if nvx
[post_mean, hpd_interval, post_var] = Extractoo(oo_, name, type);
catch
draws = getalldraws(ip);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, 1, options_.mh_conf_sig);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, options_.mh_conf_sig);
k = estim_params_.var_exo(i,1);
name = M_.exo_names{k};
oo_ = Filloo(oo_, name, type, post_mean, hpd_interval, post_median, post_var, post_deciles, density);
@ -191,7 +191,7 @@ if nvn
for i=1:nvn
if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh)
draws = getalldraws(ip);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, 1, options_.mh_conf_sig);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, options_.mh_conf_sig);
name = options_.varobs{estim_params_.nvn_observable_correspondence(i,1)};
oo_ = Filloo(oo_, name, type, post_mean, hpd_interval, post_median, post_var, post_deciles, density);
else
@ -200,7 +200,7 @@ if nvn
[post_mean,hpd_interval,post_var] = Extractoo(oo_,name,type);
catch
draws = getalldraws(ip);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws,1,options_.mh_conf_sig);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws,options_.mh_conf_sig);
name = options_.varobs{estim_params_.nvn_observable_correspondence(i,1)};
oo_ = Filloo(oo_,name,type,post_mean,hpd_interval,post_median,post_var,post_deciles,density);
end
@ -230,7 +230,7 @@ if ncx
for i=1:ncx
if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh)
draws = getalldraws(ip);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws,1,options_.mh_conf_sig);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws,options_.mh_conf_sig);
k1 = estim_params_.corrx(i,1);
k2 = estim_params_.corrx(i,2);
name = sprintf('%s,%s', M_.exo_names{k1}, M_.exo_names{k2});
@ -247,7 +247,7 @@ if ncx
[post_mean,hpd_interval,post_var] = Extractoo(oo_, NAME, type);
catch
draws = getalldraws(ip);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, 1, options_.mh_conf_sig);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, options_.mh_conf_sig);
k1 = estim_params_.corrx(i,1);
k2 = estim_params_.corrx(i,2);
name = sprintf('%s,%s', M_.exo_names{k1}, M_.exo_names{k2});
@ -281,7 +281,7 @@ if ncn
for i=1:ncn
if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh)
draws = getalldraws(ip);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, 1, options_.mh_conf_sig);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, options_.mh_conf_sig);
k1 = estim_params_.corrn(i,1);
k2 = estim_params_.corrn(i,2);
name = sprintf('%s,%s', M_.endo_names{k1}, M_.endo_names{k2});
@ -296,7 +296,7 @@ if ncn
[post_mean,hpd_interval,post_var] = Extractoo(oo_, NAME, type);
catch
draws = getalldraws(ip);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, 1, options_.mh_conf_sig);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, options_.mh_conf_sig);
k1 = estim_params_.corrn(i,1);
k2 = estim_params_.corrn(i,2);
name = sprintf('%s,%s', M_.endo_names{k1}, M_.endo_names{k2});

4
matlab/estimation/PosteriorIRF.m
View File

@ -283,7 +283,7 @@ for file = 1:NumberOfIRFfiles_dsge
for k = 1:size(temp.STOCK_IRF_DSGE,1)
kk = k+kdx;
[MeanIRF(kk,j,i),MedianIRF(kk,j,i),VarIRF(kk,j,i),HPDIRF(kk,:,j,i),...
DistribIRF(kk,:,j,i)] = posterior_moments(squeeze(temp.STOCK_IRF_DSGE(k,j,i,:)),0,options_.mh_conf_sig);
DistribIRF(kk,:,j,i)] = posterior_moments(squeeze(temp.STOCK_IRF_DSGE(k,j,i,:)),options_.mh_conf_sig);
end
end
end
@ -323,7 +323,7 @@ if MAX_nirfs_dsgevar
kk = k+kdx;
[MeanIRFdsgevar(kk,j,i),MedianIRFdsgevar(kk,j,i),VarIRFdsgevar(kk,j,i),...
HPDIRFdsgevar(kk,:,j,i),DistribIRFdsgevar(kk,:,j,i)] = ...
posterior_moments(squeeze(temp.STOCK_IRF_BVARDSGE(k,j,i,:)),0,options_.mh_conf_sig);
posterior_moments(squeeze(temp.STOCK_IRF_BVARDSGE(k,j,i,:)),options_.mh_conf_sig);
end
end
end

4
matlab/estimation/conditional_variance_decomposition_ME_mc_analysis.m
View File

@ -113,11 +113,11 @@ p_hpdsup = NaN(1,length(Steps));
for i=1:length(Steps)
if options_.estimation.moments_posterior_density.indicator
[pp_mean, pp_median, pp_var, hpd_interval, pp_deciles, pp_density] = ...
posterior_moments(tmp(:,i),1,mh_conf_sig);
posterior_moments(tmp(:,i),mh_conf_sig);
p_density(:,:,i) = pp_density;
else
[pp_mean, pp_median, pp_var, hpd_interval, pp_deciles] = ...
posterior_moments(tmp(:,i),0,mh_conf_sig);
posterior_moments(tmp(:,i),mh_conf_sig);
end
p_mean(i) = pp_mean;
p_median(i) = pp_median;

4
matlab/estimation/conditional_variance_decomposition_mc_analysis.m
View File

@ -104,11 +104,11 @@ p_hpdsup = NaN(1,length(Steps));
for i=1:length(Steps)
if options_.estimation.moments_posterior_density.indicator
[pp_mean, pp_median, pp_var, hpd_interval, pp_deciles, pp_density] = ...
posterior_moments(tmp(:,i),1,mh_conf_sig);
posterior_moments(tmp(:,i),mh_conf_sig);
p_density(:,:,i) = pp_density;
else
[pp_mean, pp_median, pp_var, hpd_interval, pp_deciles] = ...
posterior_moments(tmp(:,i),0,mh_conf_sig);
posterior_moments(tmp(:,i),mh_conf_sig);
end
p_mean(i) = pp_mean;
p_median(i) = pp_median;

4
matlab/estimation/correlation_mc_analysis.m
View File

@ -92,10 +92,10 @@ for file = 1:length(ListOfFiles)
end
if options_.estimation.moments_posterior_density.indicator
[p_mean, p_median, p_var, hpd_interval, p_deciles, density] = ...
posterior_moments(tmp,1,mh_conf_sig);
posterior_moments(tmp,mh_conf_sig);
else
[p_mean, p_median, p_var, hpd_interval, p_deciles] = ...
posterior_moments(tmp,0,mh_conf_sig);
posterior_moments(tmp,mh_conf_sig);
end
if isfield(oo_,[ TYPE 'TheoreticalMoments'])
temporary_structure = oo_.([TYPE, 'TheoreticalMoments']);

8
matlab/estimation/covariance_mc_analysis.m
View File

@ -110,11 +110,11 @@ end
if options_.estimation.moments_posterior_density.indicator
[p_mean, p_median, p_var, hpd_interval, p_deciles, density] = ...
posterior_moments(tmp,1,mh_conf_sig);
posterior_moments(tmp,mh_conf_sig);
oo_.([TYPE, 'TheoreticalMoments']).dsge.covariance.density.(var1).(var2) = density;
else
[p_mean, p_median, p_var, hpd_interval, p_deciles] = ...
posterior_moments(tmp,0,mh_conf_sig);
posterior_moments(tmp,mh_conf_sig);
end
oo_.([TYPE, 'TheoreticalMoments']).dsge.covariance.Mean.(var1).(var2) = p_mean;
oo_.([TYPE, 'TheoreticalMoments']).dsge.covariance.Median.(var1).(var2) = p_median;
@ -126,11 +126,11 @@ oo_.([TYPE, 'TheoreticalMoments']).dsge.covariance.deciles.(var1).(var2) = p_dec
if options_.contemporaneous_correlation
if options_.estimation.moments_posterior_density.indicator
[p_mean, p_median, p_var, hpd_interval, p_deciles, density] = ...
posterior_moments(tmp_corr_mat,1,mh_conf_sig);
posterior_moments(tmp_corr_mat,mh_conf_sig);
oo_.([TYPE, 'TheoreticalMoments']).dsge.contemporeaneous_correlation.density.(var1).(var2) = density;
else
[p_mean, p_median, p_var, hpd_interval, p_deciles] = ...
posterior_moments(tmp_corr_mat,0,mh_conf_sig);
posterior_moments(tmp_corr_mat,mh_conf_sig);
end
oo_.([TYPE, 'TheoreticalMoments']).dsge.contemporeaneous_correlation.Mean.(var1).(var2) = p_mean;
oo_.([TYPE, 'TheoreticalMoments']).dsge.contemporeaneous_correlation.Median.(var1).(var2) = p_median;

12
matlab/estimation/pm3.m
View File

@ -193,10 +193,10 @@ if strcmp(var_type,'_trend_coeff') %two dimensional arrays
for i = 1:nvar
if options_.estimation.moments_posterior_density.indicator
[Mean(1,i),Median(1,i),Var(1,i),HPD(:,1,i),Distrib(:,1,i),Density(:,:,1,i)] = ...
posterior_moments(squeeze(stock1(SelecVariables(i),:)),1,options_.mh_conf_sig,options_.estimation.moments_posterior_density);
posterior_moments(squeeze(stock1(SelecVariables(i),:)),options_.mh_conf_sig,options_.estimation.moments_posterior_density);
else
[Mean(1,i),Median(1,i),Var(1,i),HPD(:,1,i),Distrib(:,1,i)] = ...
posterior_moments(squeeze(stock1(SelecVariables(i),:)),0,options_.mh_conf_sig);
posterior_moments(squeeze(stock1(SelecVariables(i),:)),options_.mh_conf_sig);
end
end
else %three dimensional arrays
@ -204,21 +204,21 @@ else %three dimensional arrays
for j = 1:n2
if options_.estimation.moments_posterior_density.indicator
[Mean(j,i),Median(j,i),Var(j,i),HPD(:,j,i),Distrib(:,j,i),Density(:,:,j,i)] = ...
posterior_moments(squeeze(stock1(SelecVariables(i),j,:)),1,options_.mh_conf_sig,options_.estimation.moments_posterior_density);
posterior_moments(squeeze(stock1(SelecVariables(i),j,:)),options_.mh_conf_sig,options_.estimation.moments_posterior_density);
else
[Mean(j,i),Median(j,i),Var(j,i),HPD(:,j,i),Distrib(:,j,i)] = ...
posterior_moments(squeeze(stock1(SelecVariables(i),j,:)),0,options_.mh_conf_sig);
posterior_moments(squeeze(stock1(SelecVariables(i),j,:)),options_.mh_conf_sig);
end
if filter_step_ahead_indicator
if options_.estimation.moments_posterior_density.indicator
for K_step = 1:length(options_.filter_step_ahead)
[Mean_filter_step_ahead(K_step,i,j),Median_filter_step_ahead(K_step,i,j),Var_filter_step_ahead(K_step,i,j),HPD_filter_step_ahead(:,K_step,i,j),Distrib_filter_step_ahead(:,K_step,i,j),Density_filter_step_ahead(:,:,K_step,i,j) ] = ...
posterior_moments(squeeze(stock1_filter_step_ahead(SelecVariables(i),j,:,K_step)),1,options_.mh_conf_sig,options_.estimation.moments_posterior_density);
posterior_moments(squeeze(stock1_filter_step_ahead(SelecVariables(i),j,:,K_step)),options_.mh_conf_sig,options_.estimation.moments_posterior_density);
end
else
for K_step = 1:length(options_.filter_step_ahead)
[Mean_filter_step_ahead(K_step,i,j),Median_filter_step_ahead(K_step,i,j),Var_filter_step_ahead(K_step,i,j),HPD_filter_step_ahead(:,K_step,i,j),Distrib_filter_step_ahead(:,K_step,i,j)] = ...
posterior_moments(squeeze(stock1_filter_step_ahead(SelecVariables(i),j,:,K_step)),0,options_.mh_conf_sig);
posterior_moments(squeeze(stock1_filter_step_ahead(SelecVariables(i),j,:,K_step)),options_.mh_conf_sig);
end
end
end

34
matlab/estimation/posterior_moments.m
View File

@ -1,10 +1,11 @@
function [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(xx,info,mh_conf_sig,kernel_options)
function [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(xx,mh_conf_sig,kernel_options)
% [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(xx,mh_conf_sig,kernel_options)
% Computes posterior mean, median, variance, HPD interval, deciles, and density from posterior draws.
%
% INPUTS
% xx [double] Vector of posterior draws (or prior draws ;-)
% info [integer] If equal to one the posterior density is estimated.
% mh_config_sig [double] Scalar between 0 and 1 specifying the size of the HPD interval.
% xx [double] Vector of posterior draws (or prior draws ;-)
% mh_config_sig [double] Scalar between 0 and 1 specifying the size of the HPD interval.
% kernel_options [structure] options for kernel density estimate
%
%
% OUTPUTS
@ -21,7 +22,7 @@ function [post_mean, post_median, post_var, hpd_interval, post_deciles, density]
% kernel_density_estimate.m.
%
% Copyright © 2005-2017 Dynare Team
% Copyright © 2005-2023 Dynare Team
%
% This file is part of Dynare.
%
@ -38,20 +39,9 @@ function [post_mean, post_median, post_var, hpd_interval, post_deciles, density]
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
if nargin<4
number_of_grid_points = 2^9; % 2^9 = 512 !... Must be a power of two.
bandwidth = 0; % Rule of thumb optimal bandwidth parameter.
kernel_function = 'gaussian'; % Gaussian kernel for Fast Fourrier Transform approximaton.
else
number_of_grid_points = kernel_options.gridpoints;
bandwidth = kernel_options.bandwidth;
kernel_function = kernel_options.kernel_function;
end
xx = xx(:);
xx = sort(xx);
post_mean = mean(xx);
post_median = median(xx);
post_var = var(xx);
@ -84,8 +74,18 @@ if length(xx)>9
else
post_deciles=NaN(9,1);
end
density = [];
if info
if nargout == 6
if nargin<3
number_of_grid_points = 2^9; % 2^9 = 512 !... Must be a power of two.
bandwidth = 0; % Rule of thumb optimal bandwidth parameter.
kernel_function = 'gaussian'; % Gaussian kernel for Fast Fourrier Transform approximaton.
else
number_of_grid_points = kernel_options.gridpoints;
bandwidth = kernel_options.bandwidth;
kernel_function = kernel_options.kernel_function;
end
if post_var > 1e-12
optimal_bandwidth = mh_optimal_bandwidth(xx,number_of_draws,bandwidth,kernel_function);
[density(:,1),density(:,2)] = kernel_density_estimate(xx,number_of_grid_points,...

4
matlab/estimation/variance_decomposition_ME_mc_analysis.m
View File

@ -95,10 +95,10 @@ end
if options_.estimation.moments_posterior_density.indicator
[p_mean, p_median, p_var, hpd_interval, p_deciles, density] = ...
posterior_moments(tmp,1,mh_conf_sig);
posterior_moments(tmp,mh_conf_sig);
else
[p_mean, p_median, p_var, hpd_interval, p_deciles] = ...
posterior_moments(tmp,0,mh_conf_sig);
posterior_moments(tmp,mh_conf_sig);
end
oo_.([TYPE, 'TheoreticalMoments']).dsge.VarianceDecompositionME.Mean.(var).(exo) = p_mean;

4
matlab/estimation/variance_decomposition_mc_analysis.m
View File

@ -93,10 +93,10 @@ end
if options_.estimation.moments_posterior_density.indicator
[p_mean, p_median, p_var, hpd_interval, p_deciles, density] = ...
posterior_moments(tmp,1,mh_conf_sig);
posterior_moments(tmp,mh_conf_sig);
else
[p_mean, p_median, p_var, hpd_interval, p_deciles] = ...
posterior_moments(tmp,0,mh_conf_sig);
posterior_moments(tmp,mh_conf_sig);
end
oo_.([TYPE, 'TheoreticalMoments']).dsge.VarianceDecomposition.Mean.(var).(exo) = p_mean;