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

Changed posterior_moments (global variable options_ is not needed anymore, mh_conf_sig is a new input argument) + Cosmetic changes. 

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1836 ac1d8469-bf42-47a9-8791-bf33cf982152
time-shift
adjemian 2008-05-19 13:22:57 +00:00
parent f85752893e
commit de48d0b3ee
8 changed files with 1012 additions and 989 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
matlab/GetPosteriorParametersStatistics.m
View File

@ -1,6 +1,4 @@
function get_posterior_parameters_statistics()
% function get_posterior_parameters_statistics()
% This function prints and saves posterior estimates after the mcmc
% (+updates of oo_ & TeX output).
%
@ -62,7 +60,7 @@ if np
if options_.mh_replic
Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
[post_mean, post_median, post_var, hpd_interval, post_deciles, ...
density] = posterior_moments(Draws,1);
density] = posterior_moments(Draws,1,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
@ -96,7 +94,8 @@ if nvx
for i=1:nvx
if options_.mh_replic
Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws,1);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = ...
posterior_moments(Draws,1,options_.mh_conf_sig);
k = estim_params_.var_exo(i,1);
name = deblank(M_.exo_names(k,:));
oo_ = Filloo(oo_,name,type,post_mean,hpd_interval,post_median,post_var,post_deciles,density);
@ -130,7 +129,8 @@ if nvn
for i=1:nvn
if options_.mh_replic
Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws,1);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = ...
posterior_moments(Draws,1,options_.mh_conf_sig);
name = deblank(options_.varobs(estim_params_.var_endo(i,1),:));
oo_ = Filloo(oo_,name,type,post_mean,hpd_interval,post_median,post_var,post_deciles,density);
else
@ -161,7 +161,8 @@ if ncx
for i=1:ncx
if options_.mh_replic
Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws,1);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = ...
posterior_moments(Draws,1,options_.mh_conf_sig);
k1 = estim_params_.corrx(i,1);
k2 = estim_params_.corrx(i,2);
name = [deblank(M_.exo_names(k1,:)) ',' deblank(M_.exo_names(k2,:))];
@ -200,7 +201,8 @@ if ncn
for i=1:ncn
if options_.mh_replic
Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws,1);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = ...
posterior_moments(Draws,1,options_.mh_conf_sig);
k1 = estim_params_.corrn(i,1);
k2 = estim_params_.corrn(i,2);
name = [deblank(M_.endo_names(k1,:)) ',' deblank(M_.endo_names(k2,:))];

6
matlab/PosteriorIRF.m
View File

@ -1,6 +1,4 @@
function PosteriorIRF(type)
% function PosteriorIRF(type)
% Builds posterior IRFs after the MH algorithm.
%
% INPUTS
@ -313,7 +311,7 @@ for file = 1:NumberOfIRFfiles_dsge
for k = 1:size(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(STOCK_IRF_DSGE(k,j,i,:)),0);
DistribIRF(kk,:,j,i)] = posterior_moments(squeeze(STOCK_IRF_DSGE(k,j,i,:)),0,options_.mh_conf_sig);
end
end
end
@ -352,7 +350,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(STOCK_IRF_BVARDSGE(k,j,i,:)),0);
posterior_moments(squeeze(STOCK_IRF_BVARDSGE(k,j,i,:)),0,options_.mh_conf_sig);
end
end
end

19
matlab/dsge_posterior_theoretical_covariance.m
View File

@ -1,6 +1,4 @@
function dsge_posterior_theoretical_covariance()
% function dsge_posterior_theoretical_covariance()
% This function estimates the posterior density of the endogenous
% variables second order moments.
%
@ -11,7 +9,14 @@ function dsge_posterior_theoretical_covariance()
% None.
%
% SPECIAL REQUIREMENTS
% None.
% Other matlab routines distributed with Dynare: set_stationary_variables_list.m
% CheckPath.m
% selec_posterior_draws.m
% set_parameters.m
% resol.m
% th_autocovariances.m
% posterior_moments.m
%
%
% part of DYNARE, copyright Dynare Team (2007-2008)
% Gnu Public License.
@ -41,6 +46,11 @@ if ~rows(DrawsFiles)
DrawsFiles = dir([fname '_' type '_draws*']);
end
% Get the number of stationary endogenous variables.
nvar = length(ivar);
nar = options_.ar;% Saves size of the auto-correlation function.
options_.ar = 0;% Set the size of the auto-correlation function.
@ -104,7 +114,8 @@ for i=1:nvar
tmp(i1:i2) = Covariance_matrix(:,idx(i,j,nvar));
i1 = i2+1;
end
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(tmp,1);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = ...
posterior_moments(tmp,1,options_.mh_conf_sig);
name = fieldname(i,j,vartan);
eval(['oo_.PosteriorTheoreticalMoments.dsge.covariance.mean.' name ' = post_mean;']);
eval(['oo_.PosteriorTheoreticalMoments.dsge.covariance.median.' name ' = post_median;']);

13
matlab/dsge_posterior_theoretical_variance_decomposition.m
View File

@ -1,6 +1,4 @@
function dsge_posterior_theoretical_variance_decomposition()
% function dsge_posterior_theoretical_variance_decomposition()
% This function estimates the posterior distribution of the variance
% decomposition of the observed endogenous variables.
%
@ -11,7 +9,13 @@ function dsge_posterior_theoretical_variance_decomposition()
% None.
%
% SPECIAL REQUIREMENTS
% None.
% Other matlab routines distributed with Dynare: set_stationary_variables_list.m
% CheckPath.m
% selec_posterior_draws.m
% set_parameters.m
% resol.m
% th_autocovariances.m
% posterior_moments.m
%
% part of DYNARE, copyright Dynare Team (2007-2008).
% Gnu Public License.
@ -133,7 +137,8 @@ for i=1:nvar
post_deciles = NaN(9,1);
density = NaN;
else
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(tmp,1);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = ...
posterior_moments(tmp,1,options_.mh_conf_sig);
end
eval(['oo_.PosteriorTheoreticalMoments.dsge.VarianceDecomposition.mean.' name ' = post_mean;']);

79
matlab/kernel_density_estimate.m
View File

@ -1,72 +1,73 @@
function [abscissa,f] = kernel_density_estimate(data,number_of_grid_points,bandwidth,kernel_function)
% function [abscissa,f] = kernel_density_estimate(data,number_of_grid_points,bandwidth,kernel_function)
function [abscissa,f] = kernel_density_estimate(data,number_of_grid_points,number_of_draws,bandwidth,kernel_function)
% Estimates a continuous density.
%
% INPUTS
% data: data
% number_of_grid_points: number of grid points
% bandwidth: scalar equals to 0,-1 or -2. For a value different from 0,-1 or -2 the
% function will return optimal_bandwidth = bandwidth.
% kernel_function: 'gaussian','uniform','triangle','epanechnikov',
% 'quartic','triweight','cosinus'.
% data [double] Vector (number_of_draws*1) of draws.
% number_of_grid_points [integer] Scalar, number of points where the density is estimated.
% This (positive) integer must be a power of two.
% number_of_draws [integer] Scalar, number of draws.
% bandwidth [double] Real positive scalar.
% kernel_function [string] Name of the kernel function: 'gaussian', 'triweight',
% 'uniform', 'triangle', 'epanechnikov', 'quartic',
% 'triweight' and 'cosinus'
%
% OUTPUTS
% abscissa: value on the abscissa axis
% f: density
% abscissa [double] Vector (number_of_grid_points*1) of values on the abscissa axis.
% f: [double] Vector (number_of_grid_points*1) of values on the ordinate axis,
% (density estimates).
%
% SPECIAL REQUIREMENTS
% none.
%
% REFERENCES
% A kernel density estimator is used (see Silverman [1986], "Density estimation for statistics and data analysis")
% The code is adapted from Anders Holtsberg's matlab toolbox (stixbox).
%
% part of DYNARE, copyright Dynare Team (2004-2008)
% Gnu Public License.
if size(data,2) > 1 & size(data,1) == 1
data = transpose(data);
elseif size(data,2)>1 & size(data,1)>1
if min(size(data))>1
error('kernel_density_estimate:: data must be a one dimensional array.');
else
data = data(:);
end
test = log(number_of_grid_points)/log(2);
if (abs(test-round(test)) > 10^(-12))
if (abs(test-round(test)) > 1e-12)
error('kernel_density_estimate:: The number of grid points must be a power of 2.');
end
n = size(data,1);
%% KERNEL SPECIFICATION...
%% Kernel specification.
if strcmpi(kernel_function,'gaussian')
k = inline('inv(sqrt(2*pi))*exp(-0.5*x.^2)');
kernel = inline('inv(sqrt(2*pi))*exp(-0.5*x.^2)');
elseif strcmpi(kernel_function,'uniform')
k = inline('0.5*(abs(x) <= 1)');
kernel = inline('0.5*(abs(x) <= 1)');
elseif strcmpi(kernel_function,'triangle')
k = inline('(1-abs(x)).*(abs(x) <= 1)');
kernel = inline('(1-abs(x)).*(abs(x) <= 1)');
elseif strcmpi(kernel_function,'epanechnikov')
k = inline('0.75*(1-x.^2).*(abs(x) <= 1)');
kernel = inline('0.75*(1-x.^2).*(abs(x) <= 1)');
elseif strcmpi(kernel_function,'quartic')
k = inline('0.9375*((1-x.^2).^2).*(abs(x) <= 1)');
kernel = inline('0.9375*((1-x.^2).^2).*(abs(x) <= 1)');
elseif strcmpi(kernel_function,'triweight')
k = inline('1.09375*((1-x.^2).^3).*(abs(x) <= 1)');
kernel = inline('1.09375*((1-x.^2).^3).*(abs(x) <= 1)');
elseif strcmpi(kernel_function,'cosinus')
k = inline('(pi/4)*cos((pi/2)*x).*(abs(x) <= 1)');
kernel = inline('(pi/4)*cos((pi/2)*x).*(abs(x) <= 1)');
end
%% COMPUTE DENSITY ESTIMATE... Gaussian kernel should be used (FFT).
a = min(data) - (max(data)-min(data))/3;
b = max(data) + (max(data)-min(data))/3;
abscissa = linspace(a,b,number_of_grid_points)';
d = abscissa(2)-abscissa(1);
%% Non parametric estimation (Gaussian kernel should be used (FFT)).
lower_bound = min(data) - (max(data)-min(data))/3;
upper_bound = max(data) + (max(data)-min(data))/3;
abscissa = linspace(lower_bound,upper_bound,number_of_grid_points)';
inc = abscissa(2)-abscissa(1);
xi = zeros(number_of_grid_points,1);
xa = (data-a)/(b-a)*number_of_grid_points;
for i = 1:n;
xa = (data-lower_bound)/(upper_bound-lower_bound)*number_of_grid_points;
for i = 1:number_of_draws
indx = floor(xa(i));
temp = xa(i)-indx;
xi(indx+[1 2]) = xi(indx+[1 2]) + [1-temp,temp]';
end;
xk = [-number_of_grid_points:number_of_grid_points-1]'*d;
kk = k(xk/bandwidth);
kk = kk / (sum(kk)*d*n);
f = ifft(fft(fftshift(kk)).*fft([xi ;zeros(size(xi))]));
end
xk = [-number_of_grid_points:number_of_grid_points-1]'*inc;
kk = kernel(xk/bandwidth);
kk = kk / (sum(kk)*inc*number_of_draws);
f = ifft(fft(fftshift(kk)).*fft([ xi ; zeros(number_of_grid_points,1) ]));
f = real(f(1:number_of_grid_points));

164
matlab/mh_optimal_bandwidth.m
View File

@ -1,35 +1,44 @@
function optimal_bandwidth = mh_optimal_bandwidth(data,n,bandwidth,kernel_function)
% function optimal_bandwidth = mh_optimal_bandwidth(data,n,bandwidth,kernel_function)
function optimal_bandwidth = mh_optimal_bandwidth(data,number_of_draws,bandwidth,kernel_function)
% This function gives the optimal bandwidth parameter of a kernel estimator
% used to estimate a posterior univariate density from realisations of a
% Metropolis-Hastings algorithm.
%
% INPUTS:
% data: vector with n elements
% n: number of observations
% bandwidth: scalar equal to 0,-1 or -2. For a value different from 0,-1 or -2 the
% function will return optimal_bandwidth = bandwidth
% kernel_function: 'gaussian','uniform','triangle','epanechnikov', 'quartic','triweight','cosinus'
% data [double] Vector (number_of_draws*1) of draws.
% number_of_draws [integer] Scalar, number of draws.
% bandwidth [integer] Scalar equal to 0,-1 or -2.
% bandwidth = 0 => Silverman [1986] rule of thumb.
% bandwidth = -1 => Sheather and Jones [1991].
% bandwidth = -2 => Local bandwith parameters.
% kernel_function [string] Name of the kernel function: 'gaussian', 'triweight',
% 'uniform', 'triangle', 'epanechnikov', 'quartic',
% 'triweight' and 'cosinus'
%
% OUTPUTS:
% optimal_bandwidth: optimal window width (kernel smoothing)
% optimal_bandwidth: [double] Scalar or vector, optimal window width.
%
% SPECIAL REQUIREMENTS
% M. Skold and G.O. Roberts [2003], "Density estimation for the Metropolis-Hastings algorithm".
% Silverman [1986], "Density estimation for statistics and data analysis".
% SPECIAL REQUIREMENTS:
% none.
%
% part of DYNARE, copyright Dynare Team (2004-2007)
% REFERENCES:
% [1] M. Skold and G.O. Roberts [2003], "Density estimation for the Metropolis-Hastings algorithm".
% [2] Silverman [1986], "Density estimation for statistics and data analysis".
%
% part of DYNARE, copyright Dynare Team (2004-2008)
% Gnu Public License.
%% KERNEL SPECIFICATION...
%% Kernel specifications.
if strcmpi(kernel_function,'gaussian')
% Kernel definition
k = inline('inv(sqrt(2*pi))*exp(-0.5*x.^2)');
k2 = inline('inv(sqrt(2*pi))*(-exp(-0.5*x.^2)+(x.^2).*exp(-0.5*x.^2))'); % second derivate of the gaussian kernel
k4 = inline('inv(sqrt(2*pi))*(3*exp(-0.5*x.^2)-6*(x.^2).*exp(-0.5*x.^2)+(x.^4).*exp(-0.5*x.^2))'); % fourth derivate...
k6 = inline('inv(sqrt(2*pi))*(-15*exp(-0.5*x.^2)+45*(x.^2).*exp(-0.5*x.^2)-15*(x.^4).*exp(-0.5*x.^2)+(x.^6).*exp(-0.5*x.^2))'); % sixth derivate...
% Second derivate of the kernel function
k2 = inline('inv(sqrt(2*pi))*(-exp(-0.5*x.^2)+(x.^2).*exp(-0.5*x.^2))');
% Fourth derivate of the kernel function
k4 = inline('inv(sqrt(2*pi))*(3*exp(-0.5*x.^2)-6*(x.^2).*exp(-0.5*x.^2)+(x.^4).*exp(-0.5*x.^2))');
% Sixth derivate of the kernel function
k6 = inline(['inv(sqrt(2*pi))*(-15*exp(-0.5*x.^2)+45*(x.^2).*exp(-' ...
'0.5*x.^2)-15*(x.^4).*exp(-0.5*x.^2)+(x.^6).*exp(-0.5*x.^2))']);
mu02 = inv(2*sqrt(pi));
mu21 = 1;
elseif strcmpi(kernel_function,'uniform')
@ -68,92 +77,89 @@ else
end
%% OPTIMAL BANDWIDTH PARAMETER....
if bandwidth == 0; % Rule of thumb bandwidth parameter (Silverman [1986] corrected by
% Skold and Roberts [2003] for Metropolis-Hastings).
%% Get the Skold and Roberts' correction.
if bandwidth==0 | bandwidth==-1
correction = correction_for_repeated_draws(data,number_of_draws);
else
correction = 0;
end
%% Compute the standard deviation of the draws.
sigma = std(data);
h = 2*sigma*(sqrt(pi)*mu02/(12*(mu21^2)*n))^(1/5); % Silverman's optimal bandwidth parameter.
A = 0;
for i=1:n;
j = i;
while j<= n & data(j,1)==data(i,1);
j = j+1;
end;
A = A + 2*(j-i) - 1;
end;
A = A/n;
h = h*A^(1/5); % correction
elseif bandwidth == -1; % Adaptation of the Sheather and Jones [1991] plug-in estimation of the optimal bandwidth
% parameter for metropolis hastings algorithm.
%% Optimal bandwidth parameter.
if bandwidth == 0; % Rule of thumb bandwidth parameter (Silverman [1986].
h = 2*sigma*(sqrt(pi)*mu02/(12*(mu21^2)*number_of_draws))^(1/5);
h = h*correction^(1/5);
elseif bandwidth == -1; % Sheather and Jones [1991] plug-in estimation of the optimal bandwidth parameter.
if strcmp(kernel_function,'uniform') | ...
strcmp(kernel_function,'triangle') | ...
strcmp(kernel_function,'epanechnikov') | ...
strcmp(kernel_function,'quartic');
error('I can''t compute the optimal bandwidth with this kernel... Try the gaussian, triweight or cosinus kernels.');
end;
sigma = std(data);
A = 0;
for i=1:n;
j = i;
while j<= n & data(j,1)==data(i,1);
j = j+1;
end;
A = A + 2*(j-i) - 1;
end;
A = A/n;
strcmp(kernel_function,'quartic')
error(['I can''t compute the optimal bandwidth with this kernel...' ...
'Try the gaussian, triweight or cosinus kernels.']);
end
Itilda4 = 8*7*6*5/(((2*sigma)^9)*sqrt(pi));
g3 = abs(2*A*k6(0)/(mu21*Itilda4*n))^(1/9);
g3 = abs(2*correction*k6(0)/(mu21*Itilda4*number_of_draws))^(1/9);
Ihat3 = 0;
for i=1:n;
for i=1:number_of_draws
Ihat3 = Ihat3 + sum(k6((data(i,1)-data)/g3));
end;
Ihat3 = -Ihat3/((n^2)*g3^7);
g2 = abs(2*A*k4(0)/(mu21*Ihat3*n))^(1/7);
end
Ihat3 = - Ihat3/((number_of_draws^2)*g3^7);
g2 = abs(2*correction*k4(0)/(mu21*Ihat3*number_of_draws))^(1/7);
Ihat2 = 0;
for i=1:n;
for i=1:number_of_draws
Ihat2 = Ihat2 + sum(k4((data(i)-data)/g2));
end;
Ihat2 = Ihat2/((n^2)*g2^5);
h = (A*mu02/(n*Ihat2*mu21^2))^(1/5); % equation (22) in Skold and Roberts [2003] --> h_{MH}
elseif bandwidth == -2; % Bump killing... We construct local bandwith parameters in order to remove
end
Ihat2 = Ihat2/((number_of_draws^2)*g2^5);
h = (correction*mu02/(number_of_draws*Ihat2*mu21^2))^(1/5); % equation (22) in Skold and Roberts [2003].
elseif bandwidth == -2; % Bump killing... I compute local bandwith parameters in order to remove
% spurious bumps introduced by long rejecting periods.
if strcmp(kernel_function,'uniform') | ...
strcmp(kernel_function,'triangle') | ...
strcmp(kernel_function,'epanechnikov') | ...
strcmp(kernel_function,'quartic');
error('I can''t compute the optimal bandwidth with this kernel... Try the gaussian, triweight or cosinus kernels.');
end;
sigma = std(data);
A = 0;
strcmp(kernel_function,'quartic')
error(['I can''t compute the optimal bandwidth with this kernel...' ...
'Try the gaussian, triweight or cosinus kernels.']);
end
T = zeros(n,1);
for i=1:n;
for i=1:n
j = i;
while j<= n & data(j,1)==data(i,1);
while j<= n & (data(j,1)-data(i,1))<2*eps;
j = j+1;
end;
end
T(i) = (j-i);
A = A + 2*T(i) - 1;
end;
A = A/n;
correction = correction + 2*T(i) - 1;
end
correction = correction/number_of_draws;
Itilda4 = 8*7*6*5/(((2*sigma)^9)*sqrt(pi));
g3 = abs(2*A*k6(0)/(mu21*Itilda4*n))^(1/9);
g3 = abs(2*correction*k6(0)/(mu21*Itilda4*correction))^(1/9);
Ihat3 = 0;
for i=1:n;
for i=1:number_of_draws
Ihat3 = Ihat3 + sum(k6((data(i,1)-data)/g3));
end;
end
Ihat3 = -Ihat3/((n^2)*g3^7);
g2 = abs(2*A*k4(0)/(mu21*Ihat3*n))^(1/7);
g2 = abs(2*correction*k4(0)/(mu21*Ihat3*n))^(1/7);
Ihat2 = 0;
for i=1:n;
for i=1:number_of_draws;
Ihat2 = Ihat2 + sum(k4((data(i)-data)/g2));
end;
Ihat2 = Ihat2/((n^2)*g2^5);
h = ((2*T-1)*mu02/(n*Ihat2*mu21^2)).^(1/5); % Note that h is a column vector (local banwidth parameters).
elseif bandwidth > 0
h = bandwidth;
end
Ihat2 = Ihat2/((number_of_draws^2)*g2^5);
h = ((2*T-1)*mu02/(number_of_draws*Ihat2*mu21^2)).^(1/5); % h is a column vector (local banwidth parameters).
else
disp('mh_optimal_bandwidth:: ');
error('Parameter bandwidth must be a real parameter value or equal to 0,-1 or -2.');
error('Parameter bandwidth must be equal to 0, -1 or -2.');
end
optimal_bandwidth = h;
function correction = correction_for_repeated_draws(draws,n)
correction = 0;
for i=1:n
j = i;
while j<=n & ( draws(j,1) - draws(i,1) )<2*eps;
j = j+1;
end
correction = correction + 2*(j-i) - 1;
end
correction = correction/n;

3
matlab/pm3.m
View File

@ -40,7 +40,8 @@ function pm3(n1,n2,ifil,B,tit1,tit2,tit3,tit_tex,names1,names2,name3,DirectoryNa
tmp =zeros(B,1);
for i = 1:nvar
for j = 1:n2
[Mean(j,i),Median(j,i),Var(j,i),HPD(:,j,i),Distrib(:,j,i)] = posterior_moments(squeeze(stock1(SelecVariables(i),j,:)),0);
[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);
end
end
clear stock1

69
matlab/posterior_moments.m
View File

@ -1,29 +1,28 @@
function [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(xx,info)
% function [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(xx,info)
function [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(xx,info,mh_conf_sig)
% Computes posterior mean, median, variance, HPD interval, deciles, and density from posterior draws.
%
% INPUTS
% xx: vector of posterior draws
% info=1: estimates posterior density
% info=0: computes moments
% 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.
%
%
% OUTPUTS
% post_mean: posterior mean
% post_median: median
% post_var: variance
% hpd_interval: HPD interval
% post_deciles: post deciles
% density: density
% post_mean [double] Scalar, posterior mean.
% post_median [double] Scalar, posterior median.
% post_var [double] Scalar, posterior variance.
% hpd_interval [double] Vector (1*2), Highest Probability Density interval
% post_deciles [double] Vector (9*1), deciles of the posterior distribution.
% density [double] Matrix (n*2), non parametric estimate of the posterior density. First and second
% columns are respectively abscissa and ordinate coordinates.
%
% SPECIAL REQUIREMENTS
% none
% Other matlab routines distributed with Dynare: mh_optimal_bandwidth.m
% kernel_density_estimate.m.
%
% part of DYNARE, copyright Dynare Team (2005-2008)
% Gnu Public License.
global options_
xx = xx(:);
xx = sort(xx);
@ -32,33 +31,33 @@ post_mean = mean(xx);
post_median = median(xx);
post_var = var(xx);
n = length(xx);
m = round((1-options_.mh_conf_sig)*n);
k = zeros(m,1);
jj = n-m;
for ii = 1:m
k(ii) = xx(jj)-xx(ii);
number_of_draws = length(xx);
hpd_draws = round((1-mh_conf_sig)*number_of_draws);
kk = zeros(hpd_draws,1);
jj = number_of_draws-hpd_draws;
for ii = 1:hpd_draws
kk(ii) = xx(jj)-xx(ii);
jj = jj + 1;
end
[kmin,idx] = min(k);
[kmin,idx] = min(kk);
hpd_interval = [xx(idx) xx(idx)+kmin];
post_deciles = xx([round(0.1*n) ...
round(0.2*n)...
round(0.3*n)...
round(0.4*n)...
round(0.5*n)...
round(0.6*n)...
round(0.7*n)...
round(0.8*n)...
round(0.9*n)]);
post_deciles = xx([round(0.1*number_of_draws) ...
round(0.2*number_of_draws) ...
round(0.3*number_of_draws) ...
round(0.4*number_of_draws) ...
round(0.5*number_of_draws) ...
round(0.6*number_of_draws) ...
round(0.7*number_of_draws) ...
round(0.8*number_of_draws) ...
round(0.9*number_of_draws)]);
if ~info
density = [];
else
if info
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.
optimal_bandwidth = mh_optimal_bandwidth(xx,length(xx),bandwidth,kernel_function);
[density(:,1),density(:,2)] = kernel_density_estimate(xx,number_of_grid_points,optimal_bandwidth,kernel_function);
optimal_bandwidth = mh_optimal_bandwidth(xx,number_of_draws,bandwidth,kernel_function);
[density(:,1),density(:,2)] = kernel_density_estimate(xx,number_of_grid_points,...
number_of_draws,optimal_bandwidth,kernel_function);
end