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

v4.1: Code factorization and cosmetic changes, the marginal density use a recursive approach to compute the posterior mean and covariance. 

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@2083 ac1d8469-bf42-47a9-8791-bf33cf982152
time-shift
adjemian 2008-09-18 16:28:57 +00:00
parent af34f4d529
commit 6d08e65071
2 changed files with 34 additions and 51 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

48
matlab/compute_mh_covariance_matrix.m
View File

@ -1,4 +1,4 @@
function [m0,s0] = compute_mh_covariance_matrix()
function [posterior_mean,posterior_covariance,posterior_mode,posterior_kernel_at_the_mode] = compute_mh_covariance_matrix()
% function [m0,s0] = compute_mh_covariance_matrix()
% Estimation of the posterior covariance matrix and expectation.
@ -7,9 +7,10 @@ function [m0,s0] = compute_mh_covariance_matrix()
% None.
%
% OUTPUTS
% o m0 [double] (n*1) vector, posterior expectation of the parameters.
% o s0 [double] (n*n) matrix, posterior covariance of the parameters
% (computed from previous metropolis hastings).
% o posterior_mean [double] (n*1) vector, posterior expectation of the parameters.
% o posterior_covariance [double] (n*n) matrix, posterior covariance of the parameters (computed from previous metropolis hastings).
% o posterior_mode [double] (n*1) vector, posterior mode of the parameters.
% o posterior_kernel_at_the_mode [double] scalar.
%
% SPECIAL REQUIREMENTS
% None.
@ -33,6 +34,7 @@ function [m0,s0] = compute_mh_covariance_matrix()
global M_ options_ estim_params_
n = estim_params_.np + ...
estim_params_.nvn+ ...
estim_params_.ncx+ ...
@ -44,30 +46,32 @@ MhDirectoryName = CheckPath('metropolis');
load([ MhDirectoryName '/' M_.fname '_mh_history.mat'])
FirstMhFile = record.KeepedDraws.FirstMhFile;
FirstLine = record.KeepedDraws.FirstLine;
FirstLine = record.KeepedDraws.FirstLine;
TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
params = zeros(1,n);
oldlogpo2 = -Inf;
posterior_kernel_at_the_mode = -Inf;
posterior_mean = zeros(n,1);
posterior_mode = NaN(n,1);
posterior_covariance = zeros(n,n);
offset = 0;
m0 = zeros(n,1);
s0 = zeros(n,n);
for n = FirstMhFile:TotalNumberOfMhFiles
for b = 1:nblck
load([ MhDirectoryName '/' M_.fname '_mh' int2str(n) '_blck' int2str(b) '.mat'],'x2','logpo2');
[tmp,idx] = max(logpo2);
if tmp>oldlogpo2
oldlogpo2 = tmp;
params = x2(idx,:);
for b=1:nblck
first_line = FirstLine;
for n = FirstMhFile:TotalNumberOfMhFiles
%for b = 1:nblck
load([ MhDirectoryName '/' M_.fname '_mh' int2str(n) '_blck' int2str(b) '.mat'],'x2','logpo2');
[tmp,idx] = max(logpo2);
if tmp>posterior_kernel_at_the_mode
posterior_kernel_at_the_mode = tmp;
posterior_mode = x2(idx,:);
end
[posterior_mean,posterior_covariance,offset] = recursive_moments(posterior_mean,posterior_covariance,x2(first_line:end,:),offset);
first_line = 1;
end
[m0,s0,offset] = recursive_moments(m0,s0,x2(FirstLine,:),offset);
end
FirstLine = 1;
end
xparam1 = params';
hh = inv(s0);
fval = oldlogpo2;
xparam1 = posterior_mode';
hh = inv(posterior_covariance);
fval = posterior_kernel_at_the_mode;
save([M_.fname '_mh_mode.mat'],'xparam1','hh','fval');

37
matlab/marginal_density.m
View File

@ -46,42 +46,21 @@ TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8);
TODROP = floor(options_.mh_drop*TotalNumberOfMhDraws);
MU = zeros(1,npar);
SIGMA = zeros(npar,npar);
lpost_mode = -Inf;
fprintf('MH: I''m computing the posterior mean and covariance... ');
[posterior_mean,posterior_covariance,posterior_mode,posterior_kernel_at_the_mode] = compute_mh_covariance_matrix();
fprintf('MH: I''m computing the posterior mean... ');
for n = FirstMhFile:TotalNumberOfMhFiles
for b = 1:nblck
load([ MhDirectoryName '/' M_.fname '_mh' int2str(n) '_blck' int2str(b) '.mat'],'x2','logpo2');
MU = MU + sum(x2(ifil:end,:));
lpost_mode = max(lpost_mode,max(logpo2));
end
ifil = 1;
end
MU = MU/((TotalNumberOfMhDraws-TODROP)*nblck);
xparam1 = MU';
MU1 = repmat(MU,MAX_nruns,1);
%% lpost_mode is the value of the log posterior kernel at the mode.
fprintf(' Done!\n');
fprintf('MH: I''m computing the posterior covariance matrix... ');
ifil = FirstLine;
for n = FirstMhFile:TotalNumberOfMhFiles
for b = 1:nblck
load([MhDirectoryName '/' M_.fname '_mh' int2str(n) '_blck' int2str(b) '.mat'],'x2');
x = x2(ifil:end,:)-MU1(1:size(x2(ifil:end,:),1),:);
SIGMA = SIGMA + x'*x;
end
ifil = 1;
end
SIGMA = SIGMA/((TotalNumberOfMhDraws-TODROP)*nblck);%<=== Variance of the parameters (ok!)
MU = transpose(posterior_mean);
SIGMA = posterior_covariance;
lpost_mode = posterior_kernel_at_the_mode;
xparam1 = posterior_mean;
hh = inv(SIGMA);
fprintf(' Done!\n');
%% save the posterior mean and the inverse of the covariance matrix
%% (usefull if the user wants to perform some computations using
%% the posterior mean instead of the posterior mode ==> ).
save([M_.fname '_mean.mat'],'xparam1','hh','SIGMA');
%% end%Save.
disp(' ');
disp('MH: I''m computing the posterior log marginale density (modified harmonic mean)... ');
detSIGMA = det(SIGMA);